NASA CR-Z335 



INVESTIGATION OF NONLINEAR INVISCID 

AND VISCOUS FLOW EFFECTS 

IN THE ANALYSIS OF DYNAMIC STALL 

by Peter Crimi 

Prepared by 

AVCO SYSTEMS DIVISION 
Wilmington, Mass. 
for Langley Research Center 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. • FEBRUARY 1974 



1. Report No. 2. Government Accession No. 

NASA CR-2335 

3. Recipient's Catalog No. 

4. Title and Subtitle 

INVESTIGATION OF NONLINEAR INVISCID AND VISCOUS 
FLOW EFFECTS IN THE ANALYSIS OF DYNAMIC STALL 

5. Report Date 

February 19T4 

6. Performing Organization Code 

7. Author(s) 

Peter Crimi 

8. Performing Organization Report No. 

10. Work Unit No. 

9. Performing Organization Name and Address 

AVCO Systems Division 
Wilmington, MA 

11. Contract or Grant No. 

NAS 1-11245 

13. Type of Report and Period Covered 

Contractor Report 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 1 

Washington, DC 20546 

14. Sponsoring Agency Code 

15. Supplementary Notes 


Final Report 

16. Abstract 

A recently developed method for analyzing unsteady airfoil stall was refined by including 
nonlinear effects in the representation of the inviscid flow. Certain other aspects of the 
potential -flow model were reexamined and the effects of varying Reynolds number on stall 
characteristics were investigated. Refinement of the formulation improved the representation 
of the flow and chordwise pressure distribution below stall, but substantial quantitative 
differences between computed and measured results are still evident for sinusoidal pitching 
through stall. Agreement is substantially improved by assuming the growth rate of the 
dead-air region at the onset of leading -edge stall is of the order of the component of the 
free stream normal to the airfoil chordline. The method predicts the expected increase in 
the resistance to stalling with increasing Reynolds number. Results indicate that a given 
airfoil can undergo both trailing-edge and leading-edge stall under unsteady conditions. 


17. Key Wqrds (Suggested by Apthor(s) I 18. Distribution Statement 

Dynamic stall analysis, oscillating airfoil, 

Reynolds number, leading-edge stall Unlimited - Unclassified 


19. Security Oassif. (of this report) 20. Security Classif. (of this page) 21. No. of Pages 22. Price* 

Unclassified Unclassified 112 $4 - 25 


For sale by the National Technical Information Service, Springfield, Virginia 22151 






















INVESTIGATION OF NONLINEAR INVISCID 

AND VISCOUS FLOW EFFECTS IN THE 

ANALYSIS OF DYNAMIC STALL 

by Peter Crimi 
Avco Systems Division 


SUMMARY 


A recently developed method for analysing unsteady 
airfoil stall was refined by including nonlinear effects 
in the representation of the inviscid flow. Certain other 
aspects of the potential-flow model were reexamined and 
the effects of varying Reynolds number on stall character- 
istics were investigated. Refinement of the formulation 
Improved the representation of the flow and chordwise 
pressure distribution below stall, but substantial quanti- 
tative differences between computed and measured results 
are still evident for sinusoidal pitching through stall. 
Agreement is substantially Improved by assuming the growth 
rate of the dead-air region at the onset of leading-edge 
stall is of the order of the component of the free stream 
normal to the airfoil chordline. The method predicts the 
expected increase in the resistance to stalling with in- 
creasing Reynolds number. Results indicate that a given 
airfoil can undergo both trailing-edge and leading-edge 
stall under unsteady conditions . 



INTRODUCTION 


Unsteady airfoil stall is a problem of particular 
concern in the design and operation of helicopter rotors. 
Periodic stall and unstall of the blades at high advance 
ratio cause severe oscillatory control loads, increased 
vibration levels and, under some circumstances, a torsional 
aeroelastic instability, effectively limiting aircraft per- 
formance (Ref. 1). As a result, the problem has been the 
subject of considerable research (Refs. 2-5, for example). 

A method was recently developed for analyzing dynamic 
stall of a helicopter rotor blade. The method, which is 
described in Ref. 6, employs a model for each of the basic 
flow elements involved in the unsteady stall of a two- 
dimensional airfoil in incompressible flow. The interactions 
of these elements are analyzed by forward integration in 
time, with the aid of a digital computer. Results are in 
good qualitative agreement with measured loads, both dynamic 
lift overshoot and unstable moment variation being in clear 
evidence in the computed loading. A number of approximations 
were employed in the formulations which caused substantial 
quantitative differences, however. 

One of the approximations most open to question was the 
use of a linearized representation of the potential flow. 

This study was directed, first, to refining the model of the 
potential flow by introducing second-order terms. Certain 
other aspects of the model for a stalled airfoil were also 
considered. The method was then used to investigate the 
effects of varying Reynolds number and to determine the 
relative importance of leading-edge and trailing-edge stall 
under unsteady conditions. 
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airfoil semichord, m 

p p 

moment coefficient, C m Q ^ = m/(2 p u b ) 

p 

normal-force coefficient, C n = n/( p U b) 

p 

pressure coefficient, C p = 2 (p - p^ )/( pU ) 

airfoil chord, m 

reduced frequency, k = oj b/u 

lift per unit span, N/m 

length of dead-air region, m 

moment per unit span about quarter-chord, N 

force normal to chord line per unit span, N/m 

p 

static pressure, N/m 

p 

free-stream static pressure, N/m 

fluid speed at airfoil surface, m/s 

Reynolds number based on length indicated by 
subscript 

leading-edge radius, m 

airfoil section thickness distribution, m 
time, s 

free-stream speed, m/s 

fluid speed external to boundary layer 

foil-fixed coordinates, with origin at midchord 

angle of attack, deg or rad 

vortex strength, m/s 

boundary layer displacement thickness, m 



9 pitch angle/ deg or rad • * ; . 

P 

P fluid density, kg/i/ 

° source strength, m/s '?■' i : • 

p 

r wall shear, N/m 

/ perturbation velocity potential ! m /s 

w ' frequency of pitch oscillation,' rad/s 




6 



DESCRIPTION OP BASIC- METHOD 


The previously developed method for analyzing dynamic 
stall is outlined briefly below. Details can be found in 
Ref. 6. 

When the flow is attached (Figure la), the flow elements 
represented are: (1) a laminar boundary layer extending 

from the stagnation point over the leading edge; (2) a 
leading-edge separation bubble .(if separation occurs prior ... 
to transition); (3) a turbulent boundary layer from the re- 
attachment point of the leading-edge bubble (or the transi- 
tion point) to the trailing edge; and (4) a potential flow 
over the airfoil, including the effects of a vortical wake 
generated by the variation in time of the circulation about 
the airfoil. If the airfoil undergoes leading-edge stall 
(Figure lb), the flow elements modelled are: (1) a laminar 

boundary layer to the point of separation; (2) a laminar 
constant-pressure shear layer to the point of transition; 

(3) a turbulent constant-pressure shear layer; (4) a turbu- 
lent pressure-recovery region; and (5) a potential flow over 
the airfoil and external to the viscous mixing region, again 
including a vortical wake. If trailing-edge stall occurs 
(Figure lc), the flow elements represented are: (1) the 

laminar boundary layer; (2) the leading-edge bubble (if 
laminar separation occurs prior to transition); (3) the 
turbulent boundary layer; (4) a turbulent constant-pressure 
shear layer; (5) a turbulent pressure-recovery region; and 
(6) a potential flow with vortical wake. These elements 
were formulated as follows . 


Potential Flow 

Given the airfoil section characteristics and motions, 
together with the distribution of pressure in the dead-air 
region if the airfoil is stalled, the flow and pressure over 
the airfoil must be determined to compute the integrated 
load and analyze the boundary layer. The problem was formu- 
lated by Imposing linearized boundary conditions of flow 
tangency and pressure, using a perturbation velocity poten- 
tial derived from source and vortex distributions. The 
resulting coupled set of singular Integral equations is 
solved by casting the singularity distributions in series 
form and solving for the unknown coefficients by imposing 
boundary conditions at prescribed points. 
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Boundary Layer 


Because the relative importance of the individual 
elements of the boundary layer flow as they affect dynamic 
stall could not be established in advance, the representation 
in Ref. 6 was made as general as possible. The method of 
finite differences for unsteady flow, with variable step 
size in both streamwise and normal directions, was employed 
with the error in each finite-difference approximation of 
second order. The Cebeci-Smith eddy-viscosity model (Ref. 

7) is used to compute turbulent shear. 


Dead-Air Region 

The function of the model of the dead-air region is to 
define the streamwise distribution of pressure in that 
region, given the locations of the separation and recovery 
points and the pressure at the recovery point. The dead-air 
region is assumed to consist, in the most general case, of a 
laminar constant-pressure free shear layer from separation 
to transition, a turbulent constant-pressure mixing region, 
and a turbulent pressure-recovery region. The laminar shear 
layer is analyzed by the method of Ref. 8, assuming quasi- 
steady flow. The turbulent mixing and pressure-recovery 
regions are analyzed using the steady-flow momentum integral 
and first moment equations . Profile parameters in those 
regions are assumed to be universal functions of a dimension- 
less streamwise coordinate, with those functions derived 
from an exact viscous-inviscid interaction calculation. 
Matching of approximate solutions for the mixing and 
pressure-recovery regions at their Interface completes the 
analysis . 


Leading-Edge Bubble 

The leading-edge bubble on an unstalled airfoil is 
analyzed using the same basic relations employed for the 
dead-air region. Given the boundary-layer parameters at 
separation, the length of the bubble and the amount of pres- 
sure rise possible, for that length, in the pressure recovery 
region, are computed. That pressure rise is > compared with 
the rise in pressure in the potential flow over the length 
of the bubble. If the latter is greater than the former, 
the bubble is assumed to have burst, and the stall process 
is Initiated. 



Loading Calculation Procedure 

Calculations proceed by forward integration in time, 
given the blade motions as a function of time. If, at a 
given instant, the airfoil is not stalled, the potential- • in- 
flow is computed, and the boundary layer and leading-edge, 
bubble' are analyzed to check for bubble bursting. If the . 
airfoil is stalled, the pressure distribution in the dead- .. 
air region is computed, the potential flow evaluated, and 
the boundary layer is analyzed to locate the separation. • 
poiht . The last two steps are repeated iteratively until 
assumed and computed separation points agree. Rate of . • 
growth of the dead-air region -is determined from an- estimate.,.! 
of the rate of fluid entrainment derived from the potential- 
flow solution. In the case of leading-edge stall,, unstall 
is determined by, first postulating its occurrence and analyz- 
ing the leading-edge bubble which would then form to a seer.-: op 
tain whether that event did in fact occur. During unstall; -I; 
the dead-air region is washed off the airfoil . . The rate^ of;;;.(': 
wash-off is normally taken to be the free-stream speed. 

There is some indication, though, that the . rate should be ; .. 
considerably less than that value, as is "discussed subse- 
quently . 



FORMULATION OF POTENTIAL FLOW TO SECOND ORDER 

The following specifically concerns the derivation of 
the perturbation velocity potential to second order for 
unsteady, attached flow. The formulations derived are also 
used to compute flow and loading when the airfoil is 
stalled. However, in the latter case, while the solution 
is uniformly valid to first order, strictly consistent 
accounting of second-order terms was not attempted. A 
complete development to second order could not be justified 
without a corresponding refinement of the analysis of the 
dead-air region, which would be outside the scope of the 
present study. 

The problem of an airfoil in unsteady flow, including 
nonlinear effects, has been treated by several investigators. 
Representative of these studies is the work of Gieslng 
(Ref. 9)3 who used a finite-element approach to obtain 
numerical solutions for arbitrary transient motions, and 
Chen and Wirtz (Ref. 10), who employed an expansion of the 
acceleration potential to obtain integrated loading to 
second order for oscillatory pitching and plunging. 

The approach taken here was dictated primarily by the 
requirement for compatibility with the formulation of the 
linearized problem in Ref. 6. The velocity potential and 
boundary conditions were systematically expanded, following 
the general procedure used in Ref. 11 for the steady problem, 
as follows . 

The perturbation velocity potential can be written in 
the form 

b 

0 (x, y, t) = f a ($, t) tan -1 d * 

-b 

7 (*, t) m (V(x-f ) 2 + y 2 ) (1) 


where coordinates (x, y f ) are fixed to the airfoil, with 
origin at midchord. The surface of the airfoil, of chord 
2b, is located at y = ± T (x), - b < x < b (consideration 
has been limited to symmetric airfoils). Coordinate x 0 in 
Eq. (1) locates the end of the shed vortex wake. 
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Wake displacement effects have been omitted. Their 
contribution to the boundary conditions is of third order. 
While there is a second-order effect from wake displacement 
on the pressure at the airfoil surface, the contribution is 
symmetric and so does not affect integrated load. In any 
case, that term has reduced frequency as a factor, so for 
all practical purposes it can be regarded as third order. 

The potential is required to satisfy 


Voo + d 0/ 3 y 
Uoo + 9 0/ d x 


( y = - t (x) 

j -b < x < b 


( 2 ) 


where Uoo and Voo are apparent free-stream components, 
including effects of foil motions. 

The source and vortex distribution strengths are taken 
to be sums of terms of ascending order in angle of attack 
or a similar small parameter, i.e.. 


a = 


7 = 


0 2 + ••• 

7 -^ + 7 2 + 


If the derivatives of 0 in Eq. (2) are expanded in Taylor 
series about y = 0, terms of like order are assembled, and 
symmetric and antisymmetric contributions' separated, in the 
usual manner, it results that 


o 1 = 2 U T ' (x) 

! jr° (S, t) d ^ 

2T f x ~ 

-b 


w (x, t) 



\ 


\ 
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T ($) dt 

x -r 



2 u a 

* d x 


b 




*2 t) d$ 

x - $ 


- T7 (T V 


( 4 ) 


where U is free-stream speed and w is the downwash imposed 
by foil incidence and motions. 

The second of Eqs . (3) is, of course, the one solved 
in Ref. 6, using a Glauert-type trigonometric expansion of 
Tjl . Comparing Eqs. (3) and (4), it is seen that >2 can 
be computed using the same procedure as for 7^, with 
- (T 7 j ) 1 substituted for w. 

Some difficulty was encountered in implementing 
Eqs. (4), because T was originally formulated as. a rather 
slowly converging trigonometric series. As a result, the 
derivative of T 7 i, in series form, did not converge satis- 
factorily. The difficulty was resolved by replacing the 
original series for T with the following more rapidly con- 
verging one: 


T 


b sin 0 


1 (1 - cos e) —2. + sin 0 V t sin n 0 

2 ' b Lmd ,n 

n=l 


where r Q is leading-edge radius, and x = b cos 0. 


' ; • The same basic procedure as the one developed by 
Lighthill for steady flow (Ref. 11) was followed in deriving 
the formula for the flow at the surface, q. The x-component 
of q, denoted qx, is of the order of free-stream speed U, 
while the y-component, q y, is of first order in whatever 
expansion parameter one chooses to use. Thus, 


q = q. 


VT 


+ (%/<i x Y 
„ 2 


1 H v 


*X 


= q x + 2 


1 V 

U 


+ higher order terms 


( 5 ) 



while 


q x = u cos e p ± + y z ) + u s 


+ 


L 

-i. / 

2 IT J 

-b 


( ff-j_ + a 2 ) d $ 

_____ 


+ UTT + third-order terms 


x 


q y = U sin 0 p 


■ 2 ir / X - 

-b 


± v„ t UT ' 


7 V d £ 

F 3 

+ second-order terms 


wliere 9p is pitch angle and us and v s are contributions from 
the source distribution representing the dead-air region 
/when the airfoil is stalled (see Ref. 6). If these expres- 
sions for q x and q y are substituted in Eq. (5), a term which 
is singular, at the leading edge due to a factor (b + x) -1 /? 
is obtained, just as in the linear’ approximation . However, 
there is also a second-order term which is singular due to 
a factor (b + x)*” 1 , namely U LTT" + \ (T') 2 ], which is 
approximately equal to -. 25 Ur 0 /(b + x) near the leading 
edge.. A uniformly valid approximation for q is obtained, 
using Lighthill's procedure, by subtracting off the singular 
part of the offending term and multiplying by the factor r-r 


/ b + x A 2 

\b + x > r 0 /2 J • 

which restores the complete expression to its original form, 
to second order, some distance from the leading edge, and 
makes the result finite at the leading edge. The complete 
expression for q, uniformly valid to second order, with 
attached flow, is then 



q 




b + x 
b + x + r 


7?) 


U cos ft 


± ( 7 


J_ + 1 2 ) + u s + UTT 


X /• ( » 1+ » r ,) <U 

+ 2” / rn 


+ 


1 

2 U 


U sin e p 



. .. - . U r o f • - . . 

+ 4 (b + x) j . 

The pressure coefficient on the airfoil is computed from 

. Cj. (x, ± T(x), t) - i - (A) 2 - i ~Ji 4 ( x -°h fc ) (6) 

r - ’’r ■' , , ' 

The derivative of 0 in Eq. (6) is computed by the same formu- 
las as in Ref. 6, except that the second-order corrections 
are added to the source and vortex strenths . A term of 
second order which derives -from the Taylor expansion of 0 
about y = 0 has been omitted from Eq. (6) for the same 
reasons wake displacement effects were neglected. 
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RESULTS OP COMPUTATIONS 


All calculations were performed for a modified NACA 
0012 airfoil section. The section and a list of offsets 
are shown in Fig. 2. Unless otherwise noted, chordal 
Reynolds number was 2 million. 


Preliminary Calculations 

Computations were performed for transient pitching at 
a dimensionless rate 4b/U of .025 and a chordal Reynolds 
number of 3 million, in order to determine whether quasi- 
steady flow could be assumed in the analysis of the laminar 
and turbulent boundary layers, and so effect a substantial 
savings in computer storage requirements and running time. 
Results for the laminar boundary layer are shown in Pig. 3, 
in which are plotted the external flow magnitude u e , dis- 
placement thickness 5* and wall shear r v as a function of 
distance along the airfoil surface from the stagnation 
point. Results obtained by omitting time derivatives are 
seen to be virtually identical to those using the complete 
boundary- layer equations. 

It should be noted here that the point of vanishing 
wall shear is not generally coincident with the separation 
point in unsteady laminar flow. It was shown in Ref. 12, 
though, that a sufficient condition for identifying the 
separation point is that the boundary layer equations be- 
come invalid downstream of the true separation point. In 
the method of analyzing the boundary layer used here, there 
is a good Indication of when that event occurs, since the 
solution is obtained by iteration on the wall shear. When 
that iteration diverges, then presumably there Is no solu- 
tion to the boundary layer equations and unsteady separation 
has occurred. It was found that, for the results of Fig. 3, 
the points of vanishing wall shear and separation are 
nonetheless coincident, for all practical purposes, since 
the iteration on wall shear diverged at the first streamwise 
station downstream of the point where the wall shear went to 
zero . 


A slight difference between the quasi-steady and un- 
steady result was discernable near separation of the turbu- 
lent boundary layer downstream of reattachment of the 
leading-edge bubble, as shown in Pig. 4. However, the 
difference in separation points, 3.5 percent of chord, is 
less than the acceptable error of 3% used in the overall 
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Figure 2 AIRFOIL SECTION 


Iteration for trailing-edge stall. The assumption of' quasi 
steady flow appears to be justified, then, at least for 
those flow conditions and airfoil motions normally experi- 
enced by a helicopter rotor blade,. That assumption was 
therefore employed for all, subsequent calculations. 


Effect of Second-Order Terms 

In order to assess' the importance of the added terms 
in the potential flow model, a number of calculations were 
performed for direct comparison with results obtained using 
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DISTANCE FROM STAGNATION POINT, semichords 


Figure 3 COMPARISON OF RESULTS OF LAMINAR 
BOUNDARY LAYER ANALYSES 






Figure 4 COMPARISON OF RESULTS OF TURBULENT 
BOUNDARY LAYER ANALYSES 
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a linearized potential-flow formulation. Cases of both 
transient and sinusoidal pitching motions were considered. 

Time histories of normal force and moment coefficient 
and length of the dead-air region for a linearly varying 
pitch angle are shown in Pig. 5. The refined solution is 
seen to produce somewhat less dynamic overshoot of normal 
force than the linear formulation. However, the relative 
overshoot is about the same, the maximum static c n for the 
linearized case being somewhat larger than when second-order 
terms are included. 

A comparison of the variation of static normal-force 
and moment coefficients with angle of attack, generated by 
series of transient pitch calculations, is shown in Pig. 6 . 
Including second-order terms is seen to reduce both the 
maximum Cn and the stall angle of attack, and increase the 
slope of the C n curve below stall. The increase in 
dCn/da is due to nonlinear thickness effects. The theo- 
retical value of the slope for the airfoil analyzed is 6.82 
per radian, while the computed value is 6 . 76 ; the slight 
difference is believed to be due to numerical inaccuracies. 

Including second-order terms significantly improved 
the representation of the pressure distribution. The pre- 
dicted moment coefficient below stall using the nonlinear 
potential model is essentially zero, as it should be. The 
linear solution makes the center of pressure somewhat aft 
of the quarter-chord, the shift being due to the Lighthill 
correction which was used to remove the leading-edge singu- 
larity. 

Results for sinusoidal pitching at a reduced frequency 
of .13 are compared in Pig. 7, where C n and Cm c/4 are 
plotted against instantaneous pitch angle . The dashed 
curves are the static variations of the coefficients. The 
corresponding measured loading variation, from Ref. 2, is 
shown in Pig. 8 . Including second-order terms is seen to 
cause some dynamic overshoot, while none had been obtained 
with the linear model. The overshoot is still considerably 
less than what was measured, however. The nonlinear terms 
are seen to cause unstall to be initiated at a pitch angle 
somewhat above the static stall angle, while the linear 
solution has unstall starting slightly below that angle. 

Results for pitching at a reduced frequency of .26 are 
shown in Pig. 9. The comparable measured loading is shown 
in Pig. 10. The same effects of the second-order terms are 
evident at the higher frequency, but the relative differences 
are somewhat less . 
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Figure 7 COMPARISON OF COMPUTED LOADINGS DURING SINUSOIDAL 
PITCHING, k = .13 
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Figure 9 COMPARISON OF COMPUTED LOADINGS DURING SINUSOIDAL 
PITCHING, k= .26 





While a nonlinear potential-flow model improves the 
predicted loading somewhat, substantial quantitative differ- 
ences between computed and measured results are still 
evident. The analysis does serve to demonstrate that the 
errors in the original potential-flow model were not con- 
tributing to any great extent to those differences . Atten- 
tion can now be directed to other factors in improving the 
method, as discussed in the next section. 


Effect of Dead-Air Region Growth Rate 

As was noted in the previous description of the basic 
method, rate of growth of the dead-air region is computed 
from an estimate of rate of entrainment obtained from the 
potential-flow solution. That estimate, which is propor- 
tional to the net source strength in the dead-air region, 
can only be reasonably justified when the* length of the 
dead-air- region is near its steady state value.- However, a 
value for the growth rate is needed at the onset of leading- 
edge stall, when the region does not yet extend beyond the 
trailing edge, and is much less than its steady-state 
length of between two and three chords . It was beyond the 
scope of the original study to attempt a definitive analysis 
of this complex unsteady interaction of viscous..and inviscid 
flows, so the growth rate was arbitrarily made equal to the 
free stream speed at the onset of leading-edge stall . 

Subsequent calculations revealed, though, that one of 
the factors contributing to lift overshoot is an increment 
in lift induced on the aft portion of the airfoil when the 
dead-air region still terminates upstream of the trailing 
edge; This effect is discussed in more detail in Ref. 6. 

It would appear, then, that if the initial growth rate of 
the dead-air region were reduced, more overshoot would 
result . 

A series of transient pitch calculations were, carried 
out, parametrically varying initial growth rate £ 0 , to 
verify that this is the case. The loading time histories 
obtained are shown in Fig. 11. As expected^ the overshoot 
in Cn is increased substantially by reducing £ 0 . The 
time to reach steady state is, of course, increased as well. 

The greatest differences between computed and measured 
loading during oscillatory pitch are the amount of lift 
overshoot and the amounts of lift and moment hysteresis. 

It would appear, then,, that .considerable improvement in the 
agreement would result if j? 0 were reduced. The following 
arguments are put forth as justification for doing that. 
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The dead-air region consists of fluid garnered from 
the free stream .by the combined action of viscous shear 
(mostly turbulent, probably) and pressure gradients. Let 
it be assumed that this assimilation of fluid mass takes 
place primarily at the downstream end of the dead-air 
region, and that the fluid is trapped in the. funnel-shaped 
region subtending an angle 0 , as sketched below (lengths 
and velocities shown are approximate, assuming angle of 
attack a is small). 



Let f be the fraction of fluid entering the funnel which 
ultimately becomes part of the dead-air region. Then the 
rate of increase of the mass of trapped air is given, 
approximately, by 


m 


= f 


P ( 


fia tan 0 


)(Ua )_ 


But m is proportional to the rate of increase of area; 
m. w piia . Equating the two rates, it follows that 

£ « ( f tan 0 ) U a 

V 

Thus, if f tan 0 is a number of order one, or less, 2 0 
is of the order of TJ a , -rather than U. 

Convincing evidence that 2 0 must be substantially 
less than U was obtained by recomputing the loading for 
oscillatory pitching with k = .13, using a value of .25U 
for k Q . Presumably the process of unstall, in which the 
free stream recaptures the -fluid in-the dead-air region, 
is similar to that of stall onset, so the rate of wash-off 
of the dead-air region during unstall was also reduced from 
U to .25U. The resulting curves of C n and Cm c/4 vs. pitch 
angle are shown in Fig. 12. The computed variation of Cn 







is now seen to be in quite good agreement with the measured 
loading. Pig. 8. Maximum normal force is predicted almost 
exactly, and the areas'Within the loops are roughly the 
same. The moment variations still differ appreciably, the 
main discrepancy being between when the airfoil starts 
pitching down and where unstall begins . There is apparently 
considerably more suction over the aft portion of the air- 
foil in the dead-air region than what is predicted. The 
assumption of quasi-steady flow in the dead-air region is 
the most likely cause for this difference. 

The load during oscillatory pitching with] and 
wash-off rate of .25U was also computed for k = . 26 , with 
the results shown in Fig. 13. Again, the agreement with 
the measured loading, Pig. 10, is much improved, the peak 
C n computed being 1.88 while the measured maximum is 1.77. 
The inadequacy in the model of the dead-air region with 9p 
decreasing is also evident at the higher frequency. 

It might be conjectured that the wash-off rate should 
be less than the initial growth rate of the dead-air region, 
since reentrainment presumably is dominated by viscous 
shear, rather than, pressure gradients. Therefore, the 
effect of changing the ratio of K 0 to wash-off rate was 
investigated by repeating the oscillatory pitching -calcula- 
tions using a wash-off rate of . 1 U, with again being 
.25U. The results for k of .13 and .26 are shown in Pigs. 

14 and 15, respectively. The loading variation is seen to 
be quite insensitive to wash-off rate, regardless of the 
reduqed frequency or the value of 2 0 , so making it equal 
to 2 0 would appear to be a reasonable approximation. 


Effect of 

Reynolds Number on Stall Characteristics 

Stall characteristics of a given airfoil are necessar- 
ily a function of Reynolds number, because-of the role 
played by the boundary layer in the stall process. Speci- 
fically, one would expect that an airfoil which is subject 
to leading-edge stall at a certain Reynolds number would 
undergo trailing-edge stall at a, much 'higher ‘Reynolds 
number, since transition would then preclude formation of 
a leading-edge bubble. At intermediate Reynolds numbers, 
presumably either type of stall could occur. .Calculations _ 
were performed to investigate the effect of varying 
Reynolds number on dynamic and static stall characteristics 
in this intermediate range of Reynolds numbers. 

• ■ , J '■ ■ 
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The loading resulting from transient pitching through 
stall was analyzed for chordal Reynolds numbers of 3 and 
6 million, for comparison with the previous results obtained 
at a Reynolds number of 2 million. Pitch angle was again 
varied linearly with time up to a prescribed value and then 
held constant. 

Increasing Reynolds number to 3 million caused a 
marked increase in the resistance to bursting of the 
leading-edge bubble. The airfoil does not undergo leading- 
edge stall for pitch angles as high as about l6 degrees, 
but does experience trailing-edge stall between about 12 
and 16 degrees. At a steady-state pitch angle of 15.8 
degrees, separation of the turbulent boundary layer has 
progressed upstream to near the quarter-chord point. 

An interesting phenomenon is encountered during 
pitching through higher angles. As the separation point of 
the turbulent boundary layer moves up the chord, the resis- 
tance of the leading-edge bubble to bursting continuously 
decreases, even though the circulation and loading on the 
airfoil are decreasing as well. The reason for this is 
that the separated region has relatively little effect on 
the flow in the immediate vicinity of the leading edge, 
even though it reduces the loading over the rest of the 
airfoil. At a sufficiently' high incidence, the bubble 
bursts and leading-edge stall ensues . Results for a case 
in which both trailing-edge and leading-edge stall occur 
are shown in Pig. l6, where the loading and the separation 
point location x s are plotted against time for pitching up 
to l8 degrees. Note that very little C n overshoot is pre- 
dicted in this case. 

As expected, a further increase in Reynolds number to 
6 x 1()6 increases the resistance to both leading-edge and 
trailing-edge stall. When the steady-state pitch angle is 
15.8 degrees, the separation point is only about ilc from 
the trailing edge. At a slightly larger pitch angle, that 
point moves rapidly up the chord. Figure 17 shows the 
loading and x s variations with time for pitching up to 
16.3 degrees. The steady-state separation point is seen to 
be about .35c from the leading edge in this case. No C n 
overshoot at all is predicted for this case. 

Leading-edge stall^-ultimately occurs at high pitch 
angles for Re c = 6 x 10, as well, but the separation point 
of the turbulent boundary layer very nearly encroaches on 
the leading-edge bubble before the bubble bursts. The dis- 
tance separating them is only about .01c. 
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Results are summarized in Pig. 18, which shows the 
computed variation of static normal-force and moment coeffi- 
cients with angle of attack for. , the -three Reyholds numbers 
considered. A flagged symbol indicates that the airfoil 
was undergoing trailing-edge stall. Results of measurements 
for different Reynolds numbers for the section analyzed are 
not available. However, data in Ref. 13 show that a regular 
0012 section at Reynolds numbers between 3 and 6 million has 
a maximum lift coefficient of about 1.6, generated at an 
angle of attack of l6 degrees, which agrees fairly well with 
the computed values of maximum C n of 1.7 at Re c = 3 x 10° 
and 1.8 with' Reo = 6 x 10°, also occurring at"' about 16 
degrees angle of attack. 

The rapid falloff in normal force with angle of attack 
at higher Reynolds numbers is quite different from the be- 
havior of thicker airfoils undergoing trailing-edge stall, 
the falloff in the latter case being more gradual (see 
Ref. 14). The reason for the sharp drop-off is apparently 
that the pressure rise is quite steep near the leading-edge, 
but relatively flat aft of midchord. Thus, the separation 
point moves rapidly forward, once Incipient separation 
occurs, until it encounters the region of steep gradients 
near midchord (note the variation of x s in Figs. l6 and 17). 
On thicker airfoils, the pressure increase along the chord 
is more uniform, allowing the separation point to stabilize 
at points closer to the trailing edge. 

One case of sinusoidal pitching through stall was 
analyzed with R ec = 6 x 10° and k , = .13. The C n and Cm c/4 
variations with pitch. angle are shown in Pig. 19. Only 
trailing edge stall occurred during the cycle. As can be 
seen from Fig. 19* the normal force did not exceed the 
maximum static value . The moment variation exhibits a 
fairly large unstable (i.e., clockwise) loop, and the - 
moment undergoes rather large positive excursions. 

Since no direct determination of the type of stall has 
been made in the tests carried out to date, it is not clear 
whether the lack of lift overshoot in the calculations is 
symptomatic of trailing-edge stall or is an indication of 
an Inadequacy in the analytic model. If the latter is the 
case, the most likely source of the problem would seem to 
be, again, the assumption of quasi-steady flow in the 
dead-air region. 







Figure 19' COMPUTED LOADING DUR 
k = .13, Re c = 6 



CONCLUDING REMARKS 


The representation of the potential flow in a 
previously developed method for analyzing dynamic stall 
has been refined by including second-order terms. The 
effects of growth rate of the dead-air region during 
leading-edge stall and the effect of Reynolds number on 
stall characteristics were also investigated. Results can 
be summarized as follows. 

Including second-order terms improves the 
representation of the flow and chordwise pres- 
sure distribution on the airfoil below stall, 
but does not appreciably reduce the large 
quantitative differences between computed and 
measured results for sinusoidal pitching through 
stall . 

The amount of lift overshoot during dynamic 
stall is a strong function of the rate of growth 
of the dead-air region at the onset of leading- 
edge stall. Much improved agreement between 
theory and test with a smaller growth rate 
indicates that the growth rate is of the order 
of the component of the free stream normal to 
the airfoil chord, but further analytical and 
experimental study of the stall onset process 
is needed. In contrast to the strong dependence 
on growth rate, the loading is quite insensitive 
to the rate at which the dead-air region is 
washed off during unstall. 

The method predicts an increase in the 
resistance to both leading-edge and trailing- 
edge stall with increasing Reynolds number, as 
expected. Computations indicate that a given 
airfoil can undergo both traillng-edge and 
leading-edge stall under unsteady conditions 
at high Reynolds numbers, and that dynamic lift 
overshoot is much less when trailing-edge stall 
occurs . The latter effect may be due to the 
assumption of quasi-steady flow in the dead-air 
region . 
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APPENDIX 


PROGRAM LISTING 


A listing of the FORTRAN coding of the computer 
program follows . The program was written in FORTRAN IV 
for use on an IBM 360/75 computer. 
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C •. - 0001. 

C PROGRAM 13 ANALYZE UNSTEADY AIRFOIL STALL .0002 

C .0003 

DIMENSION US AV { 1 , 70} , U( 1, TO , 2 ) ,CMAT I 50,5 0) 0004 

p IMENS ION RMATI130), SCALE! 300, 2) , V( 1.00 ,2 ) , , UC MOO ,3 ) , ,Y (100 ) .0005 

X, SCALSI 300) '■ „ ' . - . 0006 

DIMENSION AL AMI 30), VZ IP* 30) ,FPRES( 100) ,CAMBR (24) ,X(300) . 0007 

DOUBLE PRECISION CM AT , RMAT . 0008 

DIMENSION XGAM( 30) , XC( 300 ), SBL( 300) ,XBSIG MOO) , ' > , 0009 

. DIMENSION ASZ( 301. AS( 30,30) ,BS130,30) ,ASHZ(30) , ASH ( 30 ,30) , BSH ( 30 , 3 00 10 

10 ),AR( 30), ARHI 30),UE( 300,3) . , . .. 0011 

DIMENSION BLAM( 30), FLAM( 10) , XFL AM ( 10) ‘ 0012 

DIMENSION XUt.30) , YU.I 30) , XL (.3 0)., Y.L (30) , VZ2 (30 ).,AR2( 30 ) 0013 

COMMON /LOADS/ CR , CR HAT, THICK ( 24)',THA T (24) , NF , ACAPI30 ,3 ) , ACAP2 130 » 0014 

13),NGAM,GAMAW( 1000 ), GAMA W2( 1000 ) .XIWUOOO) ,N WAKE , XSEP » X AT T , BC AP ( 10 0015 

10, 3), NS IG, RDBB»AA,BB,AKK, SS , S SLAM , ZS, UT ,UN , UI NF , XSI G ( 100 ) , : XSIGA(10 00 16 

10 ),XS IGBl 100),NSIGA,NSIGB,DXI > . . , 0017 

EQUIVALENCE (ASH(l) , SCALS(l)) , >•. 0018 

DATA I N , MOUT , MD, NT I ME , I SEP T, I WASH/5 ,6 ,250 , O', 0,2/ 0019 

DATA PI, TIME, RENEL , USTOP /3 . 1 4159 , 0. , 6.20E4,3.0/ . . . . 0020 

DAT A FP AST, FUP, VBUB /.1,.4,1./ 0021 

DATA VDFF /!./ 0022 

DATA Flam /l. 75, 1.75,1.724,1.527,1.354,1. ,.663, .452, ^25 0023 

14, .21/ • .0024 

DATA XFLAM /- l 00. , - 1 1 . 26,- 7. 01 ,-3 .48 ,-l . 766 ,0. , 1 .888 , 4 . 0025 

103,6.77,7.19/ 0026 

DATA NSBL.NZ.NY, NCORD , LOWER , MS TOP , MDTR , NOTBL , INDV/ 0027 

1 150,100, 70, 20, 1,6, 1,1,0/ 0028 

DATA EL SI G.FRZ, ARR, AMPLU, FREQ U, ALPH1 , A LPH2 ,HE AVE , AROT , FREQF, PH I H, 00 29 

V RY1, DRY,Y( 2I.TEST.UPRIM, ISEP/ . 1 ,. 06 ,3 . 5 ,0 . , 0030 

1 0.00,.1,0.,0.,-.5,0.,0.,1.001,.OG8,.01,.001 ,.0001 , 0/ 0031 

NAMELIST /CASE/ NSBL.NZ.NY, MA XT , NGAM, NSI G , NOFF .NCORD, LOWER, MSTOP, 0032 

1M0TR.N3TBL, INDV, EL S I G. DX I ,REB , RDBB ,FR Z ,ARR, AMPLU, FREQU » ALPH1 , ALPH2 0033 

1, HEAVE, AROT, FHEQT.PHIH, R Y 1 ,DP"Y , ' Y, TE ST , UPRI M ,.XU , XL , YU , Y L ,NF , 0034 

INVOR.SVOR.HVOR, BARG, X1V0R, EMI ,TORF ,SSPA ,CMPA ,CMPAS, USTOP, RENEL , 0035 

.1 FPAST, FUP.VBUB.HD, VDFF 0036 

NGAM= 24 00 37 

. NS I G= 24 ,0038 

N F=20 .0039 

U INF= 1 . . 0040 

. NWAKE=999 , 0041 

0X1=1, E10 0042. 

READ! IN, CASE) 0043 

M X=NS BL ♦NZ-l 0044 

, N DIMC= 50 0045 

INOV* INDV+l 0046 

RY= RY l 00 47 

’ ' 0048 

NOTE - OFFSETS ARE PUT IN AS LISTED IN THEORY OF WING SECTIONS, I.E 0049 
AS A FRACTION- 3F T0TAL C«ORO, XI 6EING MEASURED FROM THE . . 0050 

LEADING EDGE. 8E SURE NF I S AN EVEN NUMBER. 0051 

; ' ----- - : ..... • : 0052 

. WRITE! MOUT, 6) 0053 

. W R I T E ( M OUT, C AS El : 0054 

GO TO (343,342), INDV 0055 
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3*2 HR 1TEI MOUT , 25) NVOR , SVOR , HVOR »BARG . XI VQR ,EMI i»TORF ,SSPA 
HVOR* HVOR **2 


3*3 


7875 


350 

2*19 

2*20 

2*21 


2*22 

2*23 


BAR G= BARG/6.2832 
RDBC = ROBB / 2. 

CALL SECT (XU.YU.XL.YL.NOFF.NF ,R0BC , THICK, CAMBR)' 

ROBB = RDBC * 2. 

DO 7875 N = 1. NF 
CAMBRIA » » 2 * CAMBR ( N) 

THICMN) = 2. * THICK(N) 

CR=SQRT I R DBB)/ 2. 

CALL THCALCl CR, THICK , CRHAT, THAT ,NF I 
write(mout.a) 

WRITE( MOUT, 7 ) AMPLU. FREQU.ALPH1 .A LPH2 .HEAVE ,AROT .FREQF.ROBB, PEB 
WRITEC MOUT ,81 

WRITE! MOUT, 9) IN, CAMBR IN 1 .THICK (N) .THAT ( N) , N=1 , NF ) 

WRITE! MOUT'i T(JTI~ CRTC1THA1 

CALL SCALI S8L.NSBL.FRZ.ARR, ROBB i 
CALL C 0 RD X I N S 8L « N Z , R D 8 B , SBL « X, KC ) 

UTU=CR 

DO 350 N= l.NF 
UTU=UTU ♦THICK INI 


UTU=2.*UTU 

DO 2*20 M= 1, MX 

IF! XC( M ) - 1. ) 2 *20,2*1 9,2* 19 

M£ND=M- 1 

GO TO 2*21 

CONT IN J E 

MX-MEND 

XATT » -100. 

MXM1=MX -1 
UE ! MX +1 , 1 1 = 1. 

EPSLE= 2 •* rXTNZl-XTNZ- 1 FT 
EPSTE=X1MX I-XIMX-l) 

Al T C= 11 .50E4ZSQR T( R EB ) ' 
IF(MO.EQ.l) GO TO 2*23 


DO 2*22 M»I',"MX“ 
SCALEIH, 1)=0. 
SCALE! M, 2Tl®0v 
DO 2*22 N = 1. NY 
Unf.N, 11*TJT" 
U(M,N, 2)=0. 
NSIGA=MSTG 
NS IGB=MS I G 
NS IGl=NSTG*T~ 
HQTR =M3TR>1 
PTTCH = ALPH1 
1 FI INDV ♦ MOTR 
UN=SIN! PirCHT” 
UT=COS! PITCH) 
NOT BL=NOTBC*I 
XMAX= l.-ELS IG 


.LE. 2) PITCH = PITCH - ALPH2 


ANGS=PI/FLOAT!NS IG) 
CALL 

X S E P= l • 1 



0056 

0057 

0058 

0059 

0060 
0061 
0062 
0063 
006* 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 
007* 

0075 

0076 

0077 

0078 

0079 

0080 
0081 
0082 
0083 
008* 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 
009* 
0095 
00.96 

0097 

0098 

0099 

0100 
0101 
0102 
0103 
010* 

0105 

0106 

0107 

0108 

0109 

0110 
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00 2430 N= l, NS IGi 0111 

XS IGBI N ) = XS IGIN ) 0112 

2430 X S I GA ( sj )-=XSIG<N ) 0113 

00 2431 N = 1, NS IG 0114 

00 2431 NU = 1, 3 0115 

2431 BCAPI N» NU ) =0. 0116 

PINT=2./FL0AT( NCORO ) * 0117 

NCP l=NCORD + l 0118 

T HX 1= 1 . 5/ DX I 0119 

NGP 1=NGAM *1 0120 

NWM 1 = NH AX E— 1 0121 

COUNT = 0 . 0122 

00 8456 N= 1« NWAKE 0123 

G AM AH ( N )=0. *• ' 0124 

GAMAW2IN) = 0. 0125 

X IW(N)=1.+C0UNT 0126 

8456 C OU NT = COUNT +0X I 0127 

ANGLE=P I/FLOATINGAM I 0128 

COUNT = 0 . 0129 

DO 1002 M = 1 , N GP 1 0130 

PHI M= COUNT* AN GL E 0131 

XGAMIM)=C0S1 PHI»r» 0132 

OOUNT = 2 . 0133 

00 1001 N= 2» NGAM 0134 

AS I M* N )=COSl DOUNT *PH I M ) 0135 

1001 DOUNT= DOUNT *1 . 0136 

1002 COUNT = C0UNT*-1. v 0137 

CALL HASHIXGAM, NGAM, TIME, ALmti*LPH2iHEAVE ,AROT ,FREQF , PHIH,U 1 NF ,CA 0138 

IMBR.NF, VZ IP, 1, 1) 0139 

1 Ft INDV .EQ.2) CALL VWASHl BARG ,H VOR » SVOR »NVOR , XI VOR » U I NF , VZ I P » XG AM , 0140 

IN GP 1, DX I ) 0141 

DO 8458 M = l, NGP 1 0142 

CMAT ( M, 11 = 1. 0143 

RMATI M)=2.*VZ IPIMr 0144 

C MAT ( M, 2 ) =XGAM (Ml 0145 

DO 8457 N = 3,NGP1 0146 

8457 CMAT(M,N)=AS(M,N-1I , 0147 

8458 CONTINUE 0148 

CALL AL SOL I NGP 1, CMAT, RMAT, NDIMC) 0149 

DO 8459 N= 1, NGPi 0150 

AC API N, 1 ) = RM ATI N ) 0151 

8459 ACAPI N, 2)=ACAP1N,1) 0152 

CALL WASH2(VZ2,XGAM, NGAM, CR, THICK, NF,ACAP,0. ,0,0. ,2.) 0153 

00 8470 M = 1, NGP 1 0154 

R M ATI M ) = VZ2IM) , 0155 

CMATIM, 1)=1. 0156 

C MAT I M, 2 ) = XGAMIMJ 0157 

DO 8470 N = 3, NGPT 0158 

8470 CMATIM, N) = ASIM, N - 1) 0159 

CALL AL SOLI NGPT, CMAT rRMAT.TI 01 MCT > 0160 

00 8471 M = 1, NGP 1 0161 

ACAP2IM, 1 r=RMArmrr 0162 

8471 ACAP2IM, 2) = ACAP2(M,,1) 0163 

DO 2784 M=1,MX " 0164 

S I GN= 1 . 0165 
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I F C H-NZ I 2774,2775,2775 0166 

2774 SIGN— SIGN 0167 

2775 CALL 0 ECAL ( UE< M, 1 ) , XC I M ) , SI GN , 01 0168 

2784 JE(M, 2>=UE(H, 1) ' 0169 

00 1004 M=2,NGAM 0170 

100* BLAM(M)*( 1.125*XGAM(M)+. 1875*1 i.*XGAM(M) )* (l.-3.*XGAM(M) ) *ALCG( 11. 0171 

l+XGAMM ) )/( l.-XGAH(M) )))/OXl 0172 

BLAMl NGP1) *-1.125/0X1 0173 

C 0174 

C INDEXING IN TINE IS CARRIED OUT AT THIS POINT. 0175 

C 0176 

9999 CONTINUE 0177 

GO TO ( 196, 195I.M0TR 0178 

195 READ! IX, 21 ALPH1, ALPH2.HE AVE 0179 

C 0180 

C NOTE - FOR READ-IN OF FOIL MOT IONS , RAKF ALPH1 » ALPHA, 0181 

C ALPH2 * ALPHA-DOT, AND HEAVE = H-DOT. 0182 

C , 0183 

196 GOTO (198,197), INDV 0184 

197 CALL ELP!T(ALPH1,ALPH2»EMI« TORF , SSPA, UINF ,0X1 ,C MPA, CM PAS) 0185 

X 1 VOR=X IVOR ♦OX I*U INF 0186 

198 N ITS*l 0187 

' T IME=T IME+DX I 0188 

PITCH * ALPH1 0189 

I F ( INDv ♦ MOTR .LE. 2) PITCH * PITCH - ALPH2*CQS(FRE QF*TI ME) 0190 

NT IME=NT IME*1 0191 

NWAK E*NT IME+2 0192 

I F ( NW A< E-99B r‘"Z0Z» 207T201~ 0193 

201 NWAKE=998 0194 

202 IF! MAX T -NT THE) 8989, 8800, 8800 0195 

8800 S AV EU* J INF 0196 

U INF=1.>AHPLU*SINCTREQU*TIRE) 0197 

UBUB*UINF*VBUB 0198 

UN«=UINF*S INIPITCHT 0199 

UT»UlNF*COSI P ITCH) 0200 

U DOT* F REQU *AM PL U*COS( FREQ UPTIME J 0201 

STEPX* • 5*DX I*(UINF*SAVEU) 0202 

DO 1003 J*2,NWAXE 0 203 

JC*NWAKE-J*2 0204 

G AH AW I J C ) x GAM AM TJC - II 0205 

G AHAM2I JC ) * GAMAH2IJC-1) 0206 

1003 X lW(JC)»XrWfJC-Tl4STEPX 0207 

IF(ISEP) 2009,2009,2007 0 208 

2007 DO 2008 N*l,ffSTG 0209 

BCAPIN, 3)*BCAP(N, 2) 0210 

2008 BCAPTNi27 « BCAP T N , 1 ) 0211 

‘ DO 4433 N*1,NSIG1 0212 

XS IGB( X )=XSTGA(N I " 0213 

4433 XSIGAIN)-XSIG(N) 0214 

PRECS*PREC * . 0215 

GO TO 2010 0216 

2009 DEADEST. — 0217 

ELDOT*JBUB 0218 

2010 DO 10U M*T,MX 0219 

UE1H,3)*UE(H, 2) 0220 
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Cl Cl Cl C Cl 


1014 UECM, 2)=UE(M, 1) 0221 

. 0EAD1= DEADL ' 0222 

- ELD1=EL00T 0223 

• AL AMC 1 1 = { 1. 12 5 + . 75*AL0G< STEPX*. 5)) /OXI 0224 

■ 00 1005 M = 2,NGP1 ' ‘ 0225 

100 5- A LAM ( M I = BL AM ( M ) ♦.75*( 1 . + ( 1 ,-XGA M 1 M I) / STE P X) *ALOG ( 11 . STE PX-XG AMI M ) 0 226 

l)/< l.-XGAMCM) )>/DXI 0227 

DO 2006 M= 1, NGP 1 ' 0228 

A CAP2 ( M , 3 ) = AC AP 2( M, 2 ) 0229 

AC AP2 ( M * 2 ) = ACAP2(M,1) 0230 

AC AP( M, 3)=ACAP(M,2) 0231 

2006 AC API M, 2 ) = AC AP ( M, 1 ) 0232 

. A FACT = 8.*C ACAP ( 1, 2)*.5*ACAP 12,2)1-2.* (AC AP 1 1 ,3 ) ♦. 5* AC AP (2 .3)) 0233 

AFAC2 =8.*(ACAP2( 1, 2 ) ♦ . 5* AC AP 2 t 2 , 2 ) 1 -2 .* CACAP2 ( l ,3 ) +. 5*AC AP2 l 2 , 3 > ) 0234 

ALPHS = VZ IPC1) • - 0235 

.CALL WASH( XGAM.NGAM, TIME, ALPH 1 , ALPH2 , HEAVE, AROT ,FREQF ,PHIH ,UI NF ,CA 0236 

• 1MBR.NF, VZ 1P,M0TR, INDV ) 0237 

IFl INDV .EQ.2I CALL VWA SH C BARG ,H VOR , SVOR ,NVOR , XI VOR, UI NF , VZ l P , XGAM ,' 0238 

INGPl, D< I ) 0239 

DO 1006 M= 1, NGP 1 0 240 

ASZ CM )= l.+2.*ALAM( M) 0241 

ASCM, l )=XGAMCM)4ALAH(M) 0242 

SUM=0. 0243 

' SUM 2=0. 0244 

NWMl = NdAKE-l 0245 

DO 4343 J =2, NWM 1 0246 

JP=J+1 0247 

OLX = XlaU)-XGAMCH) - - 0248 

- C0X=DLX/(X1W<JP)-XIWIJ)) 0249 

ALX = AL3G( tXIWC JP)-XGAM(M) l/DLX) 0250 

SUM=S JM ♦ ( GAM AW I J J GAM AW 1 J )— G AMAWC JP ) ) *COX) *ALX * 0251 

4343 SUM2=S JM2*IG AMAM2t J-) » «GA«Atr2(\n^GAMAH2UP) ) *COX) *ALX 0 252 

ELX=l;-XGAM(MI l - ■ 0253 

IFCM.EQ.l) EtX=r.“ 0254 

ALX=( 1 . -X GAM ( M ) )*ALOG( l. + STEP X/EL X > /STEP X 0255 

AR2(M) = ALAH{M)*AFAC2/3.*CSOM2-GAMAW2rZ7*ALX) /PI ■ 0256 

1006 ARC M ) =2 . *VZ IP C M ) *AL AM ( M )*AFAC T/ 3. «■ I SUM-G AMA M C2 ) *ALX) / PI 0 2 57 

0258 

THE FOLLOWING CALCULATIONS, THROUGH STATEMENT 4444, ARE PERFORMED 0259 

ONLY IF THE ATRFCrn. rS "STSTtFOT - ^ THE ATRFOTL“I S DESIGNATED TO BE 0260 

STALLED IF INTEGER ISEP IS NONZERO. 0261 

0262 

IFl ISEP) 3247,4444,3247 0263 

3 247 GO TO ( 3344, 33457, I WA^H 0 264 

3344 XSEP=XSEP*DX l*VOFF ! 0265 

I F( XS EP-XMAX1 ~3248»~734T» 3347' 0266 

3347 IWASH=2 0267 

IS EP=0 0268 

' ‘ X SEP= l . I 0269 

• ' 00 3015 K-li3‘ 0270 

DO 3015 N= 1, NS IG 0271 

3015 BCAPtN,K )=tr; 0272 

GO TO 4444 0273 

3345 IFC IN DT ) 3348,3348,3248 0274 

3348 IF(NITS-l) 3248,3349,3248 0275 
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3349 IF! IND7 . EQ . 2 ) GO TO 6349 

IFIVZIPI D-ALPHS) 6349,^348,6348 

6348 N ITS= 2 

GO TO 3248 

6349 CALL UNPOPI NGP 1, A* , ALA M, A FACT ,R MAT ,CMAT,XGAM,AS,HX,NZ,XC,UE,NDIHC, 
1 AFAC2, AR2) 

GO TO 2785 

3248 XATT=XSEP*DEADI*.5*IELD1*ELD0T)*DXI 

IF< INDr .EQ.l. ANO.XATT.LT. 1.2) XATT-1.2 
DEA0L=XATT-XSEP 
0 IFF* 1 ,-XATT 

CALL SETSXINSIG1* XSEP , XATT,XSIG,ANGS) 

DO 4434 N = l» N S IG 

4434 X BS IG t N ).* . 5*1 XSIGI N J *XS1GIN *!.,))• 

DO 30 8 6 M = 1, NGP 1 

00 3086 N = 1,NSIG • • . • 

3086 BSIM,N)=0. 

IF! DlFF-l.E-6) 5005,5006, 5006 

5005 PREC=0. 

GO TO 5007 

5006 CALL ATTPRIPREC»ASZ»AS*AR,CHAT»RHAT , NGP1 ,AR2 , NO I NC , XG AH ,UTU .PRECS ) 

5007 CALL M TXE*f FPR£5fPTCEC, mNFrUOOT, TH I CfC.NF ,X8S IG , NS IG , I NOT , DELI ,THET 
11 *REB,USEP,X4«CP 1, ROBB I 

DO 3487 M*1*NGP 1 
I F( XGAMI M J-XSEP ) 3088,3088,3089 
3089 I F ( X AT r -X GAH IT 3 1 87, 3087, 3091 
3091 CONTINUE 

ACOSXiARCOST'r T^XGAWI H) - X5EP-XATT) /DEADL) 

00 30 N= 2, NSIG 

30 BSI M, N ) = SIN f FLO ATI N- 1 J* 'ACO 5X ) 

BS ( H, 1 ) =SQRT ( ( XGAMI M )-XSEP) /I XATT-XGAMIMH ) 

3088 IF( 0IFF-1.E-61 3087,3098,3098 

3098 BSIM, 1)*BSIM, 1 ) *0 IFF** I- 1 .5 )* SORT I0EADL1 * 1 2. *01 FF* l SORT I < 1. -XGAMI H 
1 ))/ (X Arr-ltOWTHTn^. m4.* XGAHTHT-1 .-3. *XATT ) ) 

GO TO 3087 

3187 BS|M, l J*DIFF**l-l,5T*$QRTt0EA0U*I3.* XATT:-4» *XGAM (H) > 

3087 CONTINUE • > , • ■ , . 

3487 CONTINUE 

C • 

C SET-UP OF THE 5TC0N0 TFT OFT QLJ A TT 0 N S STARTS HERE. 

C 

00 4350 K* l.NSTG 

1 F( XBS I Gl K )- l . ) 4348,4349,4349 
4348 COSK*X BS 1 Gt1C~) 

S INK=SQRTt l.-COSK*COSK ) 

THETK=SRCT(CffSKI 

T ANT-S INI . 5*THETX ) /CO SI • 5*THE TK I 

ASHZlK)*TAN7*CCr9A*TIT*C0S1C)Ft I.-3.*CU SK) /UINF*THXI * (PI -THETK*SI NK* 
1C0NA*! l.+COSK )*SINK**2)/UINF 
AS H I K , 1 ) = . 5* TA SHZT TC I- TSHTJ ♦ ST NIC 
COUN T* 1 . 

C OU NT* COUNT* 1. 

4355 ASHlK,NI*SINrcaONT*THETKT*.T5*rSTNnCDWrr*i; J*THETK)/ IC0UNT*1.)-S1 
IN I I COUNT- l . )*THETK I /! COUNT- 1. ) ) /IOXI *UI NF) 


0276 

0277 

0278 

0279 

0280 
0281 
0282 

0283 

0284 

0285 

0286 

0287 

0288 

0289 

0290 

0291 

0292 

0293 

0294 

0295 

0296 

0297 

0298 

0299 

0300 

0301 

0302 

0303 

0304 

0305 

0306 
030 7 

0308 

0309 

0310 

0311 

0312 

0313 

0314 

0315 
0 316 

0317 

0318 

0319 

0320 

0321 

0322 

0323 

0324 

0325 

0326 

0327 

0328 

0329 

0330 
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noon 


GO TO 4350 

4349 AS HZ ( K ) =0 . 

DO 4359 N=1,NGAM 
4359 ASH(K,N)=0. 

4350 CONT IN JE 
CPOT= CP 1 

DP 4800 K = 1» NS IG 
CORD=X B S I G C K ) 

BSHtK, 1)=-1. 

A COS X =ARCOS I I 2 .♦CORD - X SEP- XA TT) /DEADL) 

DO 4808 N= 2, NSIG 

4808 8SH«K,N>= COS I FLOAT CN-1 ) * ACOSX) 

ARH( < ) = FPRESIK) 

IF( CORO-1.) 5008,4799,4799 

500 8 CALL EGAM It 2.NGAM, ACAP ,BCAP{ 1 ,2 ) , XSJGAll) ,XSIGA INS1GAH ) .GAMAW12) . 
1C0RD, V AL 1 ) 

CALL EGAM I( 3,NGAM,ACAP,8CAP (1 ,3) .XSIGBll) ,XSIGB INSIGB+1 ) ,GAMAW(3J , 
1CORO.VAL2) 

ARHIK >=ARHlK >♦< 2.*VAL 1-. 5*VAL 2 ) / IDX1 * UI NF ) *• 062 5*AF ACT*PI * ( 1 . *CCRD 
l) *< l.-3.*C0RD*THXI*( l.-CORD*CORD))/{DXt*UINF) 

4799 CONT IN JE 

4800 CONT I N J E 
4444 CONT IN JE 

CALCULATIONS FROM THIS POINT ON COMBINE THE 
CASES OF STALLED AND ON STALLED AIRFOILS. 

DO 6500 M* l, NGP 1 
R MAT I M ) » ARI M I 
CMATIM, l)=ASZrM) 

DO 6485 N = l,NGAH 

6485 CMATIM, N + 1T=AS«H,N) 

IF(ISEP) 6486,6500,6486 

6486 00 6499 N*t7MSTG 
NGG=N +NGP 1 

6499 CMAT I M, NGG )=BS( M,N ) 

6500 CONT IN J E 

IF(ISEP) 6502,6501,6502 

6501 NT0T=NGP1 
GO TO 6751 

6502 DO 6750 K=l,NSlG 

KK*K+NGP1 .. • 

RMATI KX )= ARHlK ) / ' 

CMATI 1>=ASHZT*T 
DO 6748 N*l, NGAM 

6748 CMATlK<,"N^l)= AS Ht*- , Ni " ' 

DO 6750 N=1,NSIG 
NGG=N*NGPT 

6750 CMATIKK,NGG)=BSHIK,N) 

NTOT=NS IG*ttGM: 

6751 CALL ALSOLINTOT, CMAT, RMAT, NOIMC) 

DO 6800 N = l, N GPT 

6800 ACAPI N, 1 ) =RM AT IN) 

I FI IS EP ) 6805*“68ZD» "6805 
6805 00 6810 N=1,NSIG 


0331 

0332 

0333 

0334 

0335 

0336 

0337 

0338 

0339 

0340 

0341 

0342 

0343 

0344 

0345 
0 346 

0347 

0348 

0349 
03 50 
0351 
0 352 

0353 

0354 

0355 

0356 

0357 

0358 

0359 

0360 

0361 

0362 

0363 
0 364 

0365 

0366 

0367 

0368 

0369 
0 370 
0 371 

0372 

0373 

0374 

0375 

0376 

0377 
0 378 

0379 

0380 

0381 

0382 
0 383 

0384 

0385 
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NGG»N*NGP1 0386 

6810 BCAPlN, 1)=RHAT(NGG) 0387 

6820 CONT IN J 6 0388 

GAM AMI l ) = GAM II ACAP , DX I ,P I ) 0389 

DDEL=UTU*GAM AMI 11 0390 

DELPUOOEL/PI 0391 

DO 6910 M * It NGP 1 0392 

CMATI M, 1) * ASZIMJ 0393 

DO 6910 N * 2. NGP 1 0396 

6910 CHAT ( M» N) = ASfM, N-I) 0395 

CALL WASH2I VZ2,XGAM,NGAM,CR, THICK, NF , ACAP ,G A MAW 1 1 ) ,ISEP,DELPI , 0396 

1 XIWI2H 0397 

DO 339 M*1,NGP1 0398 

339 RHATl M) * VZ 2( M )*AR2( M ) 0399 

CALL AIS0L(NGP1,CNAT,RHAT,NDIMC ) 0600 

00 340 M= 1» NGP I 0601 

ACAP2I H» 1 1*RM AT ( M ) 0402 

340 CONTINJE 0403 

GAMAH2( 1 ) * GAM 1| ACAP 2, DXI, PIJ ♦ODE L 0404 

DO 1785 M* 1, MX 0405 

S 1 GN= 1 • 0406 

1 FI M-NZ'L IT80»~1TB~5»TT8'5 0407 

1780 SIGN--SIGN 0408 

1785 CALL QECALIUEIM, li.XCIMr, SIGN ,1 SEPT 0609 

2785 DO 8886 1*1,2 0410 

US2=UEI 1,1) 0411 

DO 8886 N* 1, MXM 1 0412 

USl*UECMVri 0413 

JEIM, 11*1 US1^US2*UEIM^1, 1)1/3. 0414 

8886 US2*US1 1 0415 

IF{ ISEPT.EQ.l) CALL TAYLORI XSEP ,UE,XC ,X,NZ,MX) 0416 

GO TO (8351»83531,I17ASfr 0417 

8351 DO 8352 M-l.MX 0418 

8352 SCACSr*7*07 0419 

GO TO 1786 0420 

8353 CALL YSETrRYl,Yr2T,NY,Yl ~ 0421 

R Y=RY 1 0422 

DO 8354 H*I,HJT 0423 

8354 SCALSM»*0. 0424 

IF! INO/.Ea72TSNDuNTIME.LT.10) WTOT rT86 0425 

I FI INDY .EQ. 2) GO T3 8370 0426 

IF! lSEP.EQ.OrANDTVZlini}.LT.Ai:PHSj GO TO 1786 0427 

8370 CALL ST AGI MX, NY, HSTOP.MST ,DXI ,R Y,DR Y,X, Y ,UE ,UC , V ,USAV ,SCALS, ISEPk 0428 

1 MO I ~ - 0429 

LAMQ-i 0430 

XSEPS*OCSTEF 0431 

DXX*DXI 0432 

IFC ISEP.EQ.I7AfnT^rS^PT.E(T.0.*N17;NITS.EQ.Il DXX?1.E30 0433 

8367 CALL 8LC( X , Y, MST,NeND,NY»RY»DRY»DXX»REB , UPRI M»F LAM* XF LAM, TEST ,U,SC 0434 
lALt,UE,UC,V,XSEP,USEP 70 rSP,THrrriXOMER,lAMQ,MSEP,XC,USAV,SCAtS,NIT 0435 

1S.NTIME, MOt 0436 

I FIXS EP— XHA X 1 77 36, 7 735, 77 3 5 0437 

7735 IFIISEP) 1786,1786,7736 0438 

7736 DEL 1=0T3P 0439 

THETl*THETA 0440 
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INDT= 1-LAMQ 0441 

I FI .NOT .( INDT.EQ.l.AND.NOTBl.EQ.2) > GO TO 378 0442 

X S EP=XS EP S 0443 

I NOT = 0 0444 

GO TO 1786 0445 

378 WRITE(M0UT,23) X S l GI1) ,CPOT , XSEP 0446 

IH INDr ) 8462, 8462, 8463 0447 

8462 IFIISEP) 8562,8562,8563 0448 

8563 IHNITS-1) 8562,8562,8662 0449 

8662 IF(ISEPT) 7742,7742,8562 0450 

8562 CALL BJ BBI DEL 1, THET 1 ,R EB , XSEP . USE P . XC 5 ,DCP ,DE L5 , X , XC , MX ,NZ ,X5 »U5 ,U 0451 

IE, ALT C, REN EL, USTOP I 1 0452 

IFIUSEP .GT.UST0P1 USEP =USEP ♦. 0751 14* I USEP-USTOP) , 0453 

PDIFF=( USEP-U5 )*IUSEP*U5) 0454 

WR ITEl *10UT, 22 > PDI FF , DCP 0455 

IF( OCP-PDIFF) 8263,8366,8366 0456 

8263 I S EPT =0 0457 

/GO TO 8463 0458 

8366 I F ( IS EP ) 8368, 8368, 8369 0459 

8369 I F ( IS EP T J 8467,8467,8368 0460 

8467 IWASH=l 0461 

N ITS 3 2 0462 

GO TO 3344 0463 

8368 GOTO ( 8168, 17861 .NOTBL 0464 

8168 CALL REATTIUC, V, X , Y ,MX ,N Y, R Y, OR Y , UE , X5 ,DE L5 , MST , REB > 0465 

LAM 0=0 0466 

GO TO 8367 0467 

8463 1 FI IS EP ) 7741,7741,7742 " 0468 

7741 ISEP=1 0469 

N IT5=N ITS+1 0470 

I F ( INOT ) 7743,7743,7643 0471 

7643 I SEPT = l 0472 

CALL CPC ( CP 1, X SEP, l.,0) 0473 

GO TO 3248 0474 

7742 CALL EL DERI BCAP , XS IG, N SIG ,UBUB , ELDOT , SI G SUM , YMX) 0475 

IFI ISEP.EO.l.AND.ISEPT.EQ .O.AND.NITS.EQ.l) GO TO 92 1 0 0 4 76 

IFIXSEP + .5J 7841,7842, 7842 0477 

7841 EPS=EPSLE 0478 

GO TO 7843 0479 

9307 I FI XSEP-XMAX1' ~9212r921-2r783 5- 0480 

9212 CALL CPC( CP l, XSEP, 1., I SEP I 0481 

GO TO 3248 0482 

7743 IF(NITS-l) 7737,7737, 3248 0483 

7737 NITS=NITS*1 , 0484 

ELDOT = ELD l 0485 

7842 EPS=EPGT"E ” 0486 

7843 DXSEP=ABSIXSEP-XSEPS) 0487 

IFI OXSEP-EPS ) 7834*7634, 9310 0488 

7834 IFIXSEP-XMAXJ 1786,1786,7835 0489 

7835 I S E P= 0 - 0490 

ISEPT=0 0491 

00 7836 0492 

DO 7836 N=l»NSIG 0493 

7836 BCAPt N, K I =0 » 0494 

GO TO 1786 0495 
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9210 

NITS=NITS*1 



0496 


I FI N ITS .EQ.2.AND.INDT . EQ . 0) XSEPS=XSEP 



0497 


l FI N IT S -4 ) 9211,9211, 1786 



0498 

9211 

I F( XS EP-XSEPS ) 9305,9305, 93 06 



0499 

9305 

X S EP=! l.-FUP l*XSEPS+FUP*XSEP 



0500 


GO TO 9307 



0501 

9306 

XSMAX=XSEPS+FPAST*( l.-XSEPS) 



0502 


IF(XSEP.GT.XSMAX) XSEP=XSMAX 



0503 


GO TO 3248 



0504 

1786 

WRITEMOUT,20) ntihe 



0505 


X SEP=XS IGI 1 ) 



0506 


WRITEI10UT.26) X1V0R 



0507 


PITC = PITCH * 180. / PI 



0508 

209 

WRITE! MOUT, 10) T 1ME , U INF , XSEP , XATT ,P I TC 



0509 


WRITE! WOUT ,11) 



0510 


WRITE! WOUT, 12 ) ( N, XGAMIN ) , V ZI P! N ) , Aft ( N) ,ACAP!N,l) , XI W ( N ) ,GAKAW ( N) , 

0511 


1N=1» NGP 1 ) 



0512 


WRITE! 6,9005) 



0513 


DO 63 M = It NGP 1 



0514 

63 

WRITE! MOOT ,9006) M, XGAM I M) , V22 1M ) ,AR 2 CM) , 

ACAP21M.D , 

XIWiM) , 

0515 

6 

A GAM AW 2! M ) 



0516 


IFIISEP) 7432,7437, 7432 ~ 



0517 

7432 

WRITE! MOUT , 13 ) 



0518 


WRITE! MOOT, 171 IN.XBSTGIN ),FPRES!N) ,ARH! N) 

, BCAP (N ,1 ) ,N= 

l.NSIG) 

0519 


WRITE! 10UT, 14) EL DOT. , 



0520 


WRITE! MOOT , 18 ) XSIG! D ,CP0T,X4,CP0T,XATT,PREC 


0521 

7433 

WRITE1M0UT.15) 



0522 


XPC=-l. 



0523 


DO 7102 N=1,NCP 1 



0524 


CALL QECALIQEL,XPC,-r.,rS6P) 



0525 


CALL QECAL!QEU,XPC, _l.,ISEP ) 



0526 


CALL CPC<CPL,XPC,-l.nSEP ) 



0527 


CALL CPC! CPUtXPC, I..ISEP) 



0528 


IF! N— 1 ) 7 5 46,7545, 75 46 



0529 

7545 

CPL=CPJ 



0530 

7546 

dlift=cpl-cpu 



0531 


WRITE! MOOT, 16) XP C, QELtCPL,QEU, CPU, D LIFT 



0 532 

7102 

XPC=XPC*PINT 



0533 


CMPAS=CMPA 



0534 


CALL CL CMT'PTTCHt'WTOTiT SEP»THP"A , C AMORT 



0535 


I FI MD .EQ • 1 ) GO TO 9999 



0536 


DO 7950 M*I7H)T 



0537 


SCALE! H,2)=SCALEIM, 1) 



0538 


SCALE!*!, I) = SCtt:S!M7 ~ 



0539 


DO 7950 N=l» NY 



0 540 


U!M,N,2)=U(M,N,n — 



0541 

7950 

UIM.N.l )=USAV(M,N> 



0542 


GO TO 9999 



0543 

8989 

CONTINUE 



0544 


STOP 



0545 

2 

FORMAT! 3F10. 4) 



0 546 

4 

F0RHATT1HT771 



0547 

5 

FORMAT! 6F10. 4) 



0548 

6 

FORMAT! IHT750X734HSNATY7IS OF“"UNSTEADY“ArrRFOI L STALL///) 

0549 

7 

FORMAT! 8X.6HUBAR =E 13.5/7X, 7HUFREQ =E 13. 5//3 X,l IHALPHA 

ONE =£13.5/ 

0550 
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13X, l 1HALPHA TWO =E13.5/8X,6HHBAR =6 l 3. 5/ l 1 X, 3H A =El 3. 5/8 X ,6HF REQ = 0551 

. IE13.5/A8X.6HR0/B =E13. 5//9X.5HREB =E13.5///> 0552 

8 ; iF0RMAT\<:l9X» 1HN.25X.4HC1N) ,26X,4HT!N) .25X.5HTH !N)/) ' 0553 

9 } FORMAT ( I20.3E30.5) 0554 

10 ...> FORMAT ( 5X » 3HT =E 13. 5/ 5X, 3HU =E13. 5/4X.4HXS =E13.5/4X,4HX0 =E13.5/4 0555 

IX.4HPA =E13.5//J .... . ", - 0556 

11 . FORMAT! ///4X, 1HN. L1X, 1HX.14X.5HVZIX) fl2X,5HRN!Xl .12X.4HAIN) ,21X,:3F 0557 

IX 1W.14X, 5HGAMMA/J . .. . 0558 

12 FORMAT! I5,4E17.5,8X,2E17.5) 0559 

13 -FORMAT! 1H1.8X, 1HN.20X, 1HX.2IX .5HFPIX1 . 22X,5HRH!N) ,2 1 X ,4HB I N J / I 0560 

14 FORMAT! //54X..9H L-DOT =E 13.5///51X.2 7HPRE SSURE S IN SEPARATED FLCW 0561 

. 1//55X, IHX, 19X.2HCP/I 0562 

15 FORMAT! 1H1, IlX.lHX, 16X.3HQEL. 15X,3HCPL,l5X,3HQEVj,15X,3HCPU.13X.9HC 0563 

, 1 PL - CPU/ ) 0 564 

16 FORMAT! 6E18.5) i ! 0565 

17 FORNAT! I10.4E25.51 ' . ' 0566 

18 FORMAT! 3! 40X , 2E20.5/J ) ' 0567 

20: FORMAT! IH1.50X.12HTIME STEP N0I3//1 0568 

22 FORMAT (;///40X, 26HINCRE ASE in CP REQUIRED I SE13. 5//40X.26HINCREASE 0569 

ilN CP POSSIBLE ISE13.5I r 0570 

23 ..FORMAT! ///45X.23HP0TENTIAL FLOW XS =E12.4/60x,8HCP!XS) =E12.4/ 0571 

1/45X, 23HBOUNOARY LAYER XS =EI2.47 0572 

25 FORMAT! 12X.4HNV =I2,3X,3HS *E 12. 4.3X.3HH »El2. 4.3X.3HG =E12.4,3X,4 0573 

, IHX1 =E12.4//12X,4HMI ,=E 12. 4i3Xi4HWT *E 12.4 ,3 X.4HPA *E12.4///I 0574 

26 1 FORMAT! 4X.4HX1 =E13.5) 0575 

101 FORMAT! ////10XV4HCRf *E13. 5/7X .7HCRHAT =E 13. 5) 0576 

9005 FORMAT! 1H1///T 50, "SECOND ORDER SOLUTION*// 0577 

90051 4XvlHir. rrx, tHX,t4Xr5rrV2txr .12X.5HRNIX) .12X.4HAIN) .13X.3H 0578 

90052X IW.14X.5HGAMMA/1 0579 .. 

9006 FORMAT! 15, 6E17. 51 ; : . , 0580 

ENO 0581 

'• -v • .* .’»• . 




) 
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SUBROUTINE BLC(X, Y, MST.MEND ,N Y,RY, DRY, DXI ,R£B,UPRIM, FLAM, XFLAM.TES ' 05B2 

IT ,U t S CAL E.UE.UC, V,XSEP, USER ,0 I S S- THF T S , LOWE R , LA MQ.MSE P, XC*USAV,SCA 0 583 

1LS,NITS,NT1ME, MX) r -• 0584 

; - ‘ ‘ • " "0585 

PROGRAM F3R ANALYZING LAMINAR AND TURBULENT BOUNDARY LAYERS' . 0586 

BY THE METHOD OF FINITE DIFFERENCES. IF THE INTEGER LAMQ 0587 

IS GREATER THAN ZERO. THE BOUNDARY LAYER IS LAMINAR. 0588 

0589 

DIMENSION U( NX , NY, 21." OSAV (MX, NY) 0590 

DIMENSION SC ALS! 300) «' X! 300),' YJ300) , UC (100,3), V( 100,2'),' '0591 

X UEI300, 3), XCI300) J “0592 

DIMENSION SOI 100), SEI 100),SF! 1001 , VI SC ! 1 00, 2) .GRAD! 100) 0593 

DIMENSION All 00), 81 100),C(100),D( 100) ,F (100) ‘V 0594 

DIMENSION ALPHA! 100), BETA! 100), GAMHAUOO) .DELTAC100) " ' '0595 

DIMENSION SCATEl300^rT,YARl riTJO I", "VAR 2 nOOT 0596 

DIMENSION FLAM! 10), X FLAM! 10), YB11 100) , YB2 1 100) ' ' 0597 

DIMENSION CAPG! 100),CAPH( 100) .CAPJllOO) .CAPKI100) 0598 

DOUBLE PRECISION API 100) , BP 1 1 00) , CP ( 1 00) ,DP( 100) ,FP (1 00) »UP( 100 ) -I- 0599 

10 FORMAT! 1HI,41X,3«H ANALYSIS OF LAMINAR BOUNDARY LA YER///5 IX , 12)<T 1 0600 

1ME STEP N0I3//51X, 12HITERATI0N NOI 3///3X ,1 HM ,8X ,1HX ,1 3 X ,2HXC *12X , 2 0601 

lHUE,iarr6H^PP / 0X ^9 X?5 H DE L TR,9X,3H0rSPLY9Xi5HTHETA,9X,5HSHEAR,4X, 0602 

1 1HI/ ) ’ ■ ■ ' ' 0603 

11 FORMAT! !Hl,4IX,56HANXtYSIS 0F~ TURBULENT BOUNDAR Y LAYER///51X, 12HT I 0604 

1ME STEP N0I3//51X,12HITERATI0N NOl 3///3X.1HM ,8X ,1HX ,13x ,2HXC, 12X ,2 0605 

1HUE, 10X , 6H^DP/OX,'9X, 5H0ECTA «9X, 5HDI SPL'»9X,5H THE TA»9X .5HSHEAR »4X » 0606 

1 1HI/ ) 0607 

12 FORMAT! 1 4.8E1 4747733 ' 0608 

20 FORMAT I 1H1.2X, 3HM -I4//2X.3HX -E14.5//2X, 4HUE -E14.5 ,10X ,17H-I 1/ R 0609 

1 H 0 )I DP/DX) ■FlAi'irrOXi WREB” i El"4.5nOX,4HU* *E14.5///i 0610 

24 FORMAT! 2X , 25HPHYSICAL DELTA -E14.5,8X,12HDELTA STAR =E14. 0611 

1 5 , 8X , 7 HTHETA" ■ET47577TX, 2WTR AN SFOR MED DELTA *E14.5«8X, 12HDE 0612 

1LTA STAR "El 4. 5, BX, 7HTHETA -E14.5///) 0613 

21 FORMAT! 25XVTFIYTT9X, 1HU, 19X* 1HVV16X, 5HDU70Y»14X,6HNUE/NU/) 0614 

22 FORMAT! 10X.5E20. 5) 0615 

23 FORMAT I //30Xn7HSEPAXATT0N Ar X «t 13.5 , 6H, XC -E13.5) 0616 

25 FORMAT! ///40X, 12HWALL SHEAR -E14.5//) 0617 

30 FORMAT! //3Or,tTHTRANSTTT0N AT-X -EI4.5) 0618 

35 FORMAT! /20X, 35HSCAL E CHANGE - Y-MAX INCREASED FR0ME12.4.3H T0E1214 0619 

1/ ) 0620 

810 FORMAT! 10X.7HAT STEPI3,22H, THE HALL GRADIENT ISE12.4) 0621 

DATA M()lir'7S/~ 0622 

8 CON » 1.5/DXI 0623 

fcon * r.Anr.*i»m — 0624 

IF(MX.GT.l) GO TO 900 0625 

8CGN-0V 0626 

FC0N*0. 0627 

DXI-1.E30 0628 

900 CONTINUE 0629 

MAX rr*TT~ 0630 

MTRAN«-1 0631 

YSUB2-YT2) 0632 

MST 2 * MST - 2 0633 

MSTl-MSr-1 0634 

GO TO I 543, 550), LOWER 0635 
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543 IF(LAMQ) 544,544,545 0636 

544 WP!TE( HOOT, 1 1 ) NT [ME , NITS - 0637 

GO TO 550 0638 

5451 WRITEtWOUT, 10) NTIHE.NITS 0639 

550 CONTINJE 0640 

YTR = SQRTIREB) 0641 

UCI l, 1) * 0. . 0642 

Vf 1,1) = 0 . 0643 

NV = NT - 2 0644 

NVM1 = NV - 1 . . : 0645 

NVP1 = NV 1 ' 0 646 

CALL YDIFFINY, ALPHA, 8ETA, GAMMA, DELTA, SO, SEVSF ,C2 ,C3 ,C4 , Y). 0647 

DO 41 N= 1, NVP 1 0648 

V TSC(N, l ) = 1. . 0649 

41 1SC( N, 2) ■* 1. ... 0650 

00 42 M=«ST2,«ST1 0651 

L = MST1-M+2 0652 

DO 50 N = 1,NV . 0653 

50 GRAD!N+l) * SD(N*1I*UC!N4-2.L) +SEI N*-l) *UC IN* l ,U-tSF (N*l) *UC1N, L) 0654 

GPAO(l) a C2*UC(2,Ll*C3*UCt3,L|*C4*UC!4,L) . 0655 

MM=M~1 0656 

CALL PGR Amror.XfUE,' DXT,PRFSS,SA,SB . SC. 5RYSS ) 0657 

DO 456 N = l, NY , • 0658 

456 UC(N, 1)=UC(N,L ) , 0659 

CALL SETUP <LAMQ,M,NV,REB,X, Y,UC .PRESS .GRAD ,DE LT ,DI S P, THETA, V I SC. HT 0660 
IRAN) 0661 

42 CONTINJE 0662 

MEND! * HFND - V , 0663 

GRADS = GRADI 1 ) ' 0664 

GRAOSS = GRADI 11 0665 

-r 0666 

THE MAIN CALCULATION STARTS HERE. 0667 

C 0668 

DO 99 i=MSTV,KENirr _ 0669 

ITER=0 . 0670 

W ALLG=0. 0671 

MP1=M*1 0672 

DELTP = OELT/YTR i 5 0673 

DISPT = DISP*YTR 0674 

THETT = THETA *YTR 0675 

SHEAR = GRADID/YTR 0676 

IF<MAXlT.EQ.tO .ANIT. ERWS.tE.EPWl MAXIT*9S 0677 

GO TO I 561, 562), LOWER . 0678 

561 WRlTEH0UT,l:2h M»XOrtfXei-Hl it>em,HTPRESS,OELTP,DISP,THETA,SHEAR» * 0679 

1 MAXIT 0680 

t GO TO 225 ' 0681 

562 WRITE! ^OUT, 20 ) M , X ( M ) , UE I M, 1 ) »PRE S S,REB , UPRI M 0682 

, WRITE! *OUT, 24) DEL TP i OTSP , THE TA .DELT.DT SPT, THETT . 0683 

WRITE! M0UT.2 l ) ' 0684 

WRITE! MOUT ,22 ) ( Y( Ml. UC!N ,2 1 * V!N i t > ,GRAD IN) , VI SC < N,1 ) , N=1 , NVPl I 0685 

WRITE! MOUT, 25) SHEAR ! 0686 

225 I E! GR ADSS - GRAUS-"TiE- 6 ) 22 9,229,408 — 0687 

408 XSX=K l M-2 I •*•! X ! M- 1 )-X( H-2 ) )*GR AD SS/t GRADS S-GR AOS I 0688 

v IF<XSX-X(Mtr 4O9V409Y229 0689 

409 WFS=!XSX-X!M-1))/!X!M)-X!M-1) ) , 0690 
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GO TO 224 0691 

229 I F( GRADC l ) ) 227, 227,2000 0692 

2000 IF(MAXIT.LT.IO) GO TO 223 0693 

IF! EPMS.GT . EP W ) GO TO 223 0694 

HFS=0. 0695 

GO TO 224 - ' 0696 

227 WFS=GRAD$/< GRADS-GRADI II I 0697 

224 WFS 1= l.-HFS 0698 

XSEP»WFS1*XC(M-11+WFS*XC(M) 0699 

XBL*WFS1*X(M-1)4WFS*X<MJ 0700 

USEP=WFSl*UE!M-l, 1H-WFS*UE(H, 11 0701 

WFP*(XBL-X(M-2) l/(X!M- ll-X<M-2l 1 0702 

WFP1M.-WFP 0703 

DISS=DIS5S*WFP l*OISS*WFP 0704 

THETS = THETSS*WFPUTHETS*MFP 0705 

I F( THETS «LT .0,1 THETS*-THETS 0706 

I F( M.L T .MST 1+3 ) GO TO 223 0707 

WRITEH0UT.23I XBL.XSEP ' 0708 

IFILAHQ.EQ.0.AND.M.LT.MTRAN45 J LAMQ“1 0709 

GO TO 222 0710 

223 CONTINJE 0711 

IF(LAMffrWI,^0r.'802 • 0712 

802 CALL TRANSIUPRIM, PRESS, THETA, RE B «UCiNY, FLAM, XF LAM, LAMQ) ' 0713 

IF(LAMQ) 805,805,801 0714 

805 WRITE! ROUT , 30 1 XCM| < 0715 

MT RAN * M + l 0716 

801 CONTINJE 0717 

I FI V« l-DELTT 6Z0, 641 , 64 1 0718 

620 RY*RY*DRY 0719 

0720 

RESCALING CALCULATION STARTS HERE,, 0721 

0722 

00 632 N* 1 »NY 0723 

YB1IN) s TTNT 0724 

V AR1( N I « UC ( N, 2 I 0725 

632 VAR2INI * UCIN.Tl 0726 

CALL YSETIRY,YSUB2,NY,Y) 0727 

WPITEIN0UT,35T Y8KWY»,YCNYI 0728 

00 633 N-2.NVP1 0729 

YIN * TIN) 07 30 

CALL TERPI Y IN , YBIt VAR 1»NY ,UPA SI I 0731 

UC! N, 2 1 * UP AST 0732 

CALL T ERPI YIN, YBIt VAR 2, NY »UPA S2 I 0733 

633 UC(N, 3} « UPAS2 ' 0734 

CALL YDIFF!NY, ALPHA, BETA, GAMMA, DELTA. SO, SE,SF,C2.C3,C4,Y1 .0735 

I F I C AMO 1' 70TT,~7 00/ 701 " 0736 

700 00 635 N*2,NVP 1 0737 

VARUNl * VISCTNVU ‘ 0738 

635 VAR2!N) * VISC(N,2) 0739 

00 636 N-2.NVPT 0740 

YIN « Y ( N I 0741 

V ISCIN, 1) « U PASI 0743 

CALL T ERP (YIN , YBl, VAR2,NYP1 ,UPASZ) 0744 

636 V ISC! N, 2 ) * UPAS2 0745 
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701 00 637 N=2»NVP 1 0746 

VAR1IN) = VIN, 1) 0747 

637 VAR21NJ = V(N,2) ' 0748 

- DO 638 N= 2» N VP 1 . 0 749 

Y IN = Y IN ) 0750 

CALL TERPIYIN,YB1,VAR1,NVP1,UPAS1> ■' 0751 

VIN.1) = UPAS 1 0752 

CALL TERPIYIN.YBl, VAR 2 ,NVP 1 , UPA S2 ) 0753 

638 VIN, 2) = UPAS2 •. ' ' 0754 

641 CONT IN J E 0755 

0756 

RESCALING CALCULATION ENDS HERE. • ' ' - • ' 0757 

0758 

CALL PGRADIM,X,UE,DXI,PRESS,SA,SB,SC ,SR,SS) ' 0759 

0760 

RECURSION RELATIONS ARE SE T UP HERE. 0761 

• - 0762 


521 

522 


523 

510 

B20 


574 


575 


5 76 
88 


IFIMX.EQ.il GO TO 820 - 0763 

IF! SCAL EIM + 1, 1 )- 1 • ) 522, 522.521 0764 

IFISCALEl M+l, 21-1.1 522,522,523 0765 

LACKJ=1 - 0766 

F ACU 1= J E ( M ♦1 » 2)/t)Et M* 1 , t 1 0767 

FACU2=JEIM+1,3)/UE<M*1,1) j 0768 

GO TO 820 0769 

LACKU=2 • 0770 

DO 610 NN=1,NY • ' 0771 

VARKNNI = UIM + l.NN.l) 0772 

VAR2INNI * U I M ♦ 1 » NN , 2 1 0773 

CALL YS ET (SCALE! M ♦1, 1 ) , Y SUB 2, N Y , YB l J 0774 

CALL YS ET ( SCALE! M*l, 2 I , YSU8 2, NY , YB2I f; 0775 

DO 88 N = 2, NV ' ' • •• 0776 

CALL CAPS! ITER, N, CAP G, CAP H; CA PJ ,CAPK, SR, SS, SO, SE,SF, VISC, V,UC) 0777 

AIN) = -SF(N)*CAPG(N)-DELTA(N»*CAPH(NH-SF(N»*CAPJ(NI 0778 

B ( N) = BCQN+SA^CAPKTTJT^FfN T* r CAPG TN I-GA M7W IN)*CAPH I N) — SE IN) *CAPJINI 07 79 

C(N»=SD(N)*CAPG(N)-BETA(N)*CAPHIN>-SD<N»*CAPJ(N1 0780 

DIN) = -ALPHA! N)*CAPHLN) 0781 

IF(MX.EQ.l) GO TO 576 0782 

GO TO I 574,575), LACKU 1 0783 

U P AS 1 =F ACU1*UC ( N , l ) ’ 0784 

U PAS2=FACU2*UCLt(i't'7 ' 0785 

GO TO 576 0786 

Y IN = YIN) 0787 

CALL TERP(YIN,YBI,VAR1,NY,UPAS1) 0788 

CALL TERPI YIN# YB^tVAR2-,NY-,UPA S2I 0789 

FIN) = PRESS«-FC0N*14.*UPAS1-UPAS2I*CAPK(N)*(SR*UC(N,2)~SC*UC(N,3) )' 0790 

CONTINUE- 0791 

0792 

SOLUTION FOR V ELOCTTr PROF ILE STARTS HERE. ^ 0793 

0 794 


89 


DO 89 

N = 2 ,NV'~ 

0795 

APIN) 

= AIN) 

• 0796 

BPIN) 

= bcnt 

0797 

CP(N) 

= CIN) 

0798 

OPIN) 

= OCN) 

0799 

FPIN) 

= FIN) 

0800 


* 
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DO 77 N-2, NVM 1 0801 

CP(N) » CP<NI/BP<N) 0802 

O p ( N ) = 0P<N)/8P<N) 0803 

FP( N) = FP(N)/BP(N) 0804 

BP(N«-1) = BP(NM) - CP ( N ) *AP( N* 1 ) 0805 

CP( N + l ) a CP(N+1> - DP (N ) *AP{ N* 1) 0806 

77 FP( N+ 1 ) = FPIN+l) - FP(N)*AP< N*l) 0807 

U P ( NY I = uE( M*l» 1) 0808 

, UPlNVPl) = UP I NY ) 0809 

DP(NV) = <FP(NV)-UPINY)*tDPCNV) *■ C1MNV) D/BPINV) 0810 

DO 66 N-3.NV 0811 

NN-NV +-2-N ‘ 0812 

66 UP(NN) = FP(NN ) - DPI NN >*UP INN* 2) - CPI NNI*UPINN*1 1 0813 

DO 65 N-2, NY 0814 

65 UC(N. It ■ UP(N) ... • 0815 

I F I ITER) 843.841,843 0816 

841 DO 842 N=2*NVP 1 r 0817 

V ( N i 2 ) = V ( N» 1) 0818 

842 VISCIN»2)=VI$CIN»1) 0819 

D1SSS=D1SS 0820 

DISS=OISP ' 0821 

T HETSS=THET~S 0822 

T HETS-T HET A . , 0823 

GRADSS-GRADS 0824 

GRADS-GRADI II 0825 

843 DO 55 N-2 , NVP 1 0826 

55 VIN.ll - VI N- l« II-.54IYINI-YIN- 11 1 * t SA* I IX ( N,l ) *UCt N-l ,11 1 -SB* CUCC 0827 

lN,2H-JCCN-l,2tr*SC»niCCIIf3)4UC(I^T,3m 0828 

DO 56 N«1,NV /■ 0829 

56 GRAD(NM) = SDIN*l>*UCIN*2,l)>SE(N«-Tr*UC IN*1,1I-SFIN*1)*UCIN,1> 0830 

GRADll) = C2*UCI2. 1)*C34UC( 3,11 *C4*X 14,1) 0831 

CALL SETUPILAMQ,NP1,NV,REB,X, Y,UC, PRESS, GRAD.DELT.DISP, THETA, VISC, 0832 

1MTRANI : - 0833 


809 
8 30 
811 


119 

120 
814 

812 

44 


48 

99 


ITER- ITER3T' , 0834 

GO TO ( 830, 809), LOWER " 0835 

WR ITEI ROUT, 810) ITER.GRAOtl) 0836 

IFC ITER-9) 811,811,812 0837 

EPWS-EPW 0838 

EPW=ABSI GRADI ll-WALLG) 0839 

IF(WALL&-r.I T2IT,T207119 0840 

EPW-EPrf/WALLG 0841 

IFIEPW-TEST1 812,814, 814 0842 

l>l ALLG-GRADI 1 ) 0843 

GO TO 820 0844 

DO 44 N - 1 * NY 0845 

UC(N, 3) » TTCLNrZt 0846 

JCt N, 2) * UC( N, 1 ) 0847 

MAX IT* ITER 0848 

IFIMX.EQ.l) GO TO 99 0849 

SCALSIfm)-RY 0850 

00 48 N-l, NY 0851 

CONTINUE’ ’ 0853 

XSEP-l.l 0854 

USEP-UEIMX, 1 ) 0855 
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222 CONI' IN J E 
RETURN 
END 


0856 

0857 

0858 



v. ; 
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SUBROUTINE STAG! MX, NY , MSTOP.M ST ,0 XI ,R Y ,DR Y ,X , Y,UE ,UC, V, US AV ,S C ALS , 0859 

1 1 SEP* MD I 0860 

C PROGRAM FOR CALCULATING THE BOUNDARY LAYER PROFILE NEAR 0861 

C THE STAGNATION POINT 0862 

C 0863 

DIMENSION USAVIMD, NY), SCALSI1) 0864 

DIMENSION PHIZ! 24),PHiPl 24) ,E TAPI24) 0865 

DIMENSION XI 300),Y( 100),UE( 300,3) ,UC(100,3) ,V(100,2) 0866 

DIMENSION EFI 1007, EFPIIOOT 0867 

DATA ETAP /O . , .2, .4, .6 , .8, 1 ., 1. 2, 1.4, 1. 6 , 1 . 8 ,2. ,2 .2 ,2 .4 ,2 .6 ,2 .8 , 3 . 0868 

1. 3.2. 3.4. 3.6. 3. 8. 4.. 4.2. 4.4. 4. 6/ 0869 

DATA PHIZ /O., .0233, .0881, .1867, .3124, .4592, .622, .7967, .9798, 1.168 0870 

19, 1.362, 1.5578,1.7553, 1.9538, 2. 153, 2. 3526,2. 552 3,2. 7522 ,2. 952 1 ,3. 1 0871 

15 21,3.3521, 3.5521,3.7521,3.9521/ 08 72 

DATA PHIP /O., .2266, .4145, V5663,. 6859,. 7 7 7 9, . 84 67 ,. 89 68 ,.9323 , .95 6 0 8 7 3 

18. . 9732. .9839, .9905, . 9946,. 997, . 9984, . 9992 ,. 9996 ,.9998 ,.9999,1 . , 1 . 0874 

1.1.. 1./ 0875 

BAG*. 08 0876 

IF(ISEP) 10, 10,5 0877 

5 BAG*. 35 0878 

10 EFI 1 ) * 0. 0879 

EFPI 1 ) * 0. 0880 

DO 20 H-l.MX 0881 

IFIUE14, 1)) 20,20, 19 0882 

19 MSP * M 0883 

GO TO 21 0884 

20 CONTINJE 0885 

21 ASTAG * IUEIMSP+2, 1 )— UEI MSP +l,l))/IXt MSP+2 )— XIMSP+1 ) ) 0886 

I FI AST AG) 22,22,23 0887 

22 AST AG* l U El M SP, 1 )-UEI M SP- 1 , 1 ) ) /( XI MSP! -X1MSP— 1 ) ) 0888 

23 SQAS = SORT! AS TAG! 0889 

DELT = 2.6/SQAS 0890 

309 I FI DEL T -YTNY-3I J 3 11, 310, 317) 089 1 

310 RY*RY+QRY 0892 

CALL YSETIRY.YIZT.NY.Y) 0893 

GO TO 309 0894 

311 CONTINJE 0895 

DO 80 N*2,NY 0896 

YET « rtTU*STT*r 0897 

DO 33 NN*1, 24 0898 

ifiyet-etapinn n 408, 408,33 0899 

408 MARK = NN 0900 

GO TO 410 0901 

33 CONTINJE 0902 

EFCNT* YET-. 6479 0903 

EFPIN) - 1. 0904 

GO TO 80 0905 

410 FRACT * I YET-ETAPI MARK-1 ) )/ IE TAPI MARK ) -ETAP I MARK-1 ) I 0906 

FRAC1 * T.-FRXCT 0907 

EFIN) * PHIZ I M ARK-1 ) *FRAC 14PH IZ 1MARK)*FRACT 0908 

80 CONTINJE 0910 

Ml * MSP-HSTUP 0911 

M2 = MSP+MSTOP 0912 
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M=M 1-1 
50 M=M+1 

MST = M + 1 

00 71 N=1,NY 

; UC<N,3) = UC(N.2) 

J C( N« 2 ) = UEIM,1)*EFP(NI 
V(N,2) = V(N, l) 

V(N,l) * -SQAS*EF«N ) 

71 CCNTINJE 

1 F( MO . EQ.il GO TO 73 
SCALS ( H i- RY 

DO 72 9= 1 »NY 

72 J S AVI M, N ) =UC ( N • 2 ) 

73, ; IFIM-M21 50, 55,55 

55 IFIUEI 9, II -BAG I 50, 50, 81 
81 CCNTINJE 
RETURN 
END 


0913 

0914 
09 15 

• 0916 

0917 

0918 

0919 

0920 

0921 

0922 

0923 

0924 

0925 

0926 
09 27 
09 28 

0929 

0930 
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SUBROUTINE THCALC!CR, THICK, CRHAT, THAT, NFI 0931 

DIMENSION THICK! 24)»THATI24),ARG(61) 0932 

. DATA NS IMP, P I /60, 3. 14159/ 0933 

; S IGN= 1 • 09 34 

SUM- 0 . '> 0935 

;■ DO 5 N=1,NF 0936 

T HAT! N I =0 . ' v; , « ; 0937 

SIGN*-SIGN . • 0938 

5 SUM*SJM+S IGN*THlCK<N) ' * ! 0939 

CP HAT »C.R*I 3.*CR*2.*SUM ) ' : 0940 

. DT=PI/FLOAT!NSlMP| •' 0941 

ARGI 1 1*0. 0942 

APGINS IMP ♦l ) =0 • 0943 

i ANGLE*!) • 0944 

DO 10 N=2,NSIMP " 0945 

ANGLE* ANGLE+DT 0946 

S INT*$ IN! ANGLE ) \ 0 947 

C0ST=C3S! ANGLE ) 0948 

ST=0. 0949 

STPR=0. 0950 

COUNT *0. 09 51 

00 6 I* 1, NF 0952 

COUNT*COUNT*l. 0953 

TTI*COJNT*ANGLE 0954 

SINIT*SIN!TTM 0955 

ST*ST ♦THICK! II*SINIT 0956 

6 STPR=STPR*THrCRrn^rCirUNT*SINTrsiNIT-2.*COST*COS(TTm 0957 

ARG!N)*<CR*! 1 .-COST » *S INT*ST) ♦ICR* ll.-2.*C0ST» ♦ STPR! . 0958 

10 ARGINJ-! ARC! N) -CRHAT*! 1.-C0STT1 /SINT 0959 

0T2=DT+DT 0960 

ANGLE=0. 0961 

DO 15 N*2»NS IMP, 2 0962 

NPl-NTl 0963 

ANGLE»ANGLE*0T2 0964 

ANGA* AN GL E—DT~ 0965 

COUNT *0 • 0966 

DO 12 1*1, NF 0967 

COUNT* COUNTS I. 0968 

12 THAT! II*THATriT*4.*SRGTTn»STNn:DUNT»ANGA)f?.»ARGCNPn*SIN!CCUNT*AN 0969 

1GLEI 0970 

15 CONTINUE 0971 

FACT=DT 2/ ( 3.*P I ) 0972 

DO 20 N*I,NF 0973 

20 THAT! N l = FACT* THAT IN) 0974 

RETURN 0975 

END 0976 
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■ SUBROUTINE CLCM ( P ITCH, AROT , I SEP ,CMPA , CAMBR ) 

COMMON /LOADS/ CR , CR HA T , THI CK 1 24 ) , THAT (24) ,NF ,ACAP (30 ,3 ) ,ACAP2(30, 
1 3 ), NOAM. CAM AW ( l 000 ) , GAMA W2< 1000) , Xt W ( 1 000) , NWAKE » XSE P » XATT .BCAPI10 
10, 3), NS IG, ROBB, AA,BB,AKK, SS,SSLAM ,ZS,UT,UN,UI NF ,XSIG (100) , XS I GA ( 10 
10 ),XS IGB( 100),NSIGA,NSIGB,DXI . - ' , 

DIMENSION CAMBR (24), X( 2),T( 2) ,CPL I 2 ) ,CPU(2) ,C P 1 2 ) , T P (2 ) , A (2 ) , 

1 6(2), C(2) 

DATA NSIMP.PI /80, 3.14159/ 

CN = 0. - .. 

CT=0. 

C MZ =0 . 

T 1=0. 

T F= P I 

IFIISEP.EQ.l .AND. XATT.LT.l. ) TE =ARCOS< XA TT) 

99 DT=(T6-Tl)/FL0ATINSIMP) 

X ( 1)=C3S( Tl) . . 

CALL CPC( CPL ( 1 ), X( 1 ) , - 1. , I SEP ) 

CALL CP,C(CPU< 1 ),X( 1), 1. , I SEP ) 

CALL PR IMESCCP(l). TPC1), Tl.CR .THICK, CAMBR, NF ) 

SUMN=(CPL( 1)-CPU( 1))*SIN(T1) 

SUMT=( CP( 1 )“TP(”1) T*CPL (1 ) — tCPT l )*- TPt 1 ) )*CPU( 1 ) 

SUMM=SJMN*X< 1) 

T ( 2 )=T 1 

DO 10 N = 2, NS IMP, 2 
T I 1 )*T ( 2) +0T 
T ( 2 )=T ( 1 ) +DT 
DO 5 1*1,2 
X ( I ) = C3SI T( I) ) 

CALL CPCI CPL I I),X( I), -l., I SEP) 

CALL CPC< CPU( I ), X( I ) , 1. , I SEP ) 

CALL PR IMESICP I I ), TPC I ),TC I ) ,CR , THICK ,CAHBR , NF ) 

A ( I )= ( CPL ( I )-CPU< I ) )* S IN ( T( I ) ) 

B ( I > = I C PI I) -TP C T ) ) * CPL I 1 ) - 1 CP IT I *TP II) I *CP U 1 1 ) 

5 C( I) = A( I)*X( I) 

SUMN=SUMN+4.*A( i)*2.*A(?> 

SUMT=SJMT+4.*B( 1)+2.*B(2J 
10 SUMM=SJMM+4.*C( 1)*2.*C( 2) 

CN=CN+QT ♦(SUMN-AI2) )/6. 
c T = C T ♦ D T * cstm T-e t t*-.- 
CMZ=CMZ-DT*( SUMM-C ( 2 ) ) / 12. 

IF(TE.EQ.PI) GO TO 20 
T 1 = TE 
T E= PI 
GO TO 99 

20 COSP=COSI PTTCHT 
SINP=SIN(P ITCH) 

CL=CN*COSP-CT*SINP 
CD=CT *COSP+CN*S INP 
C MP A*CMZ + . 5*AR0T*CN 

WRITE! 6, 2) CL, CO,CN»CT,CHZ,CNPA ,A ROT 

l 4HCT =E13.5//39X, 5HCMZ =E1 3. 5/ 3 BX , 6HC MPA =E 13. 5 ,5X ,4H I A = E 13 .5 *1H 
1) ) 

RETURN 

END 
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* 

SUBROUTINE CPCICP.XC, SIGN, ISEPJ 1032 

COMMON /LOADS/ CR, CRHAT, THI CK ( 241 , THA T ( 24) ,NF , ACAP { 30 ,3 ) , ACAP2 ( 30 , 1033 

13 ),NGAM, GAMAMl 1000) *GAMAW2( 10 00) * XI WIlOOO) ,NWAKE ,XSEP,XATT .BCAPI 10 10 34 

10,3), NS 1G.RDBB, AA.SB.AKK, SS, $ SLAM ,ZS. UT ,UN, Ul NF , XSI G < 100) .XSIGAUO 1035 
10 ),XSIGB( lOO),NSIGA,NSIGB,DXt 1036 

DIMENSION GAMS! 3),ASUM( 30,3) 1037 

NGP1 = NG AM ♦ 1 1038 

DO 6 I * 1, 3 1039 

GAMS(I) = GAMAW1 II ♦ GAHAW2TT1 1040 

DO 6 J = l, NGP l 1041 

6 ASUMIj.I) = ACAPIJ.I) * ACAP21J,!) 1042 

RECIP*l./I UIN F*U INF ) 1043 

CALL QECAL I U, XC, S IGN, ISEP ) 1044 

CP = (U/UINF)**2 - 1. 1045 

CALL EG AM K I* N6AM,ASUM,BCAP1 1 • 1 ) » XStG 1 1 T ,'XS1 G INSIG+ 1 ) ,GAM SUl.XC, 1046 
1VAL1) 1047 

CALL EGAM It 2,NGAM, ASUM «BCAP tl ,2) , XS1GA 1 1 ) , XSIGA (NSI GA+1 ) ,GAP S(2), 1048 

1XC.VAL2) 1049 

CALL EGAM I ( 3, NGAM, ASUM »BCAP 1 1,3 ) , XSIGB (1 ) , XS IGB (NSI GB+1 ) , GAM SI3), 1050 

1XC.VAL3) 1051 

CP-CP+S I GN*RECIP»IT.5*VALl-2.*VA L2*V5»VAtD /DXI 1052 

20 CP=-CP 1053 

RETURN 1054 

END 1055 
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SUBROUTINE QECALI Q , X , S IGN , I SE P) 1056 

COMMON /LOADS/ CR, CRHAT, THICK ( 2 4) , THAT < 24) ,NF ,ACAP <30 ,3 ) , ACAP2 <30 , 1057 

13 >,NGAM,GAHAW< 1000), GAMAX2<1000),XTWl 1000) ,NNAKE »XSEP»XATT ,BCAP< 10 1058 

/ 10,3),NSIG,RDBB,AA,BB,AKK, SS, SSLAH»2S,UT»UN,UINF ,KSIG< 100) ,XSIGA<10 >1059 
10),XSIGB< 100),NSIGA,NSIGB,OXI ,, 4060 

DIMENSION SN< 30),CS« 30) 1061 

R DT-R DB8/2 • . ; : ' ''10 62 

USIG-O. - 1063 

X P- l.+X ; 1064 

XM= l.-X v 1065 

COST-X 1066 

T HETA=ARCOS< X I 1067 

NSUM-NGAM 1068 

IFINF.GT.NGAM) NSUM-NF 1069 

ANGLE-0. 1070 

DO 5 N* 1, NSUM 1071 

ANGLE-ANGLE+THETA 1072 

SN< N ) = S INI ANGLE) 1073 

5 CS < N) -COS < ANGLE ) 1074 

SINT-SNI1I 1075 

SAS-O. 1076 

S AC-ACAP1 1, 1) 1077 

DO 6 N* 1, NGAM 1078 

NP-N*l 1079 

SAS*SAS*< ACAPINP, l)«ACAP2<NP,2) )*SN<N) 1080 

6 SAC»SAC*ACAP<NP, 1)*CS(N) 1081 

COST2-COST+CD5T 1082 

ST-O. 1083 

ST 1*0 . 1084 

ST2-0. 1085 

ST3*0. 1086 

EN»0. 1087 

ENS-O. 1088 

DO 7 N-l.NF 1089 

EN-EN+1 • 1090 

ENS-ENS *S INT 1091 

ST ICK-SNI N )*THTCK< N1 1092 

ST3»ST3*STICK*EN*EN 1093 

ST-ST+SnCK ~ 1094 

ST 1-ST !♦< THICK INI ♦THATIN) )* <ENS* SNIN) -C0ST2*CSI Nl I 1095 

7 ST2-ST2*EN*THICKIN)*CS(N) 1096 

ST3«3.PC0ST*ST2-SINT*<2.*ST1*ST3) 1097 

ST2"SINT*ST?»Ct)ST2»ST ' 1098 

ALTP-SA$*.25*XP*< 3.*X- l. )*<GAMAM( 11+GAMAW2 ll) ) 1099 

ALTH-ATCAP1 1,11 FA CAP 21 1,1) 1100 

IFI ISEP.EO.O) GO T) 10 1101 

OEAOL-X ATT-XSEP ~ 1102 

CALL S I GF < X, US IG, XSEP, XATT, BC AP ,NS IG ,DEADL ,RDT) 1103 

TF(XATT.GE.ITT“GO TO TO' 1104 

IFIX.LT.XATTI GO TO 10 1105 

ALTP-XirTP ♦<.( l.»3.»XA TT -3.*R0 T -4. » X) »SQR T <DEA OV»)mrtX-XATT»ROT ) 1/1 1106 

11 . -X ATT ) * *‘1. 5)»BC AP < 1, 1) 1107 

10 ALT P- . 5 *S I GN* AL TP FUS IG ♦UINF*T TCRFCRHA T)* IXH-X) ♦ ST1 ♦ < 1 . -X«- X *X ) *CR** 1108 

12+< CR*<M+SINT*ST)*< 3.*X*CR+ST3) ) 1109 
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SUBROJT INE S IGF ( X , U S I G , XSEP , X AT T , BCAP ,NS IG ,DEADL , ROT ) 

DIMENSION BCAP I 100, 31 ’• 

I F ( X.LE.XSEP .OR. X.GE.XATTI GO TO AO '' ' 

COSB=( 2.*X-XATT-XSEP J/DEAOL 
BET A=ARCOS( COSB) 

ANGLE = 0. . " " ' 

SUMC=0 . 

00 10 N=2tNSIG 

ANGLE=ANGLE+BETA . ..... . , . ... . . 

10 1 SUMC=SUMC>BCAP (N» 1 )*COSlANGLE ) ‘ *''■'* . . 

JS IG= .5*1 SUB C- BCAP ( 1» 1) ) 

RETURN 

AO XBAR=(2.*X-XATT-XSEP)/0EADL 

SUMAN0=XBAR-X8AR*SQRT ( ABSI 1 .- XSAR** (-21 ) ) ‘ 

‘ F ACTOR= SUMAND 

SUM=BCAP( 2, 1 )*SUMANO 
00 45 N = 3. NS IG 

•SUMAND=SUMAND*FACTQR . 

45 SUM=SJM+BCAP(N; 1)*$UMAND v; / ' 
i US IG= SUM*. 5*BCAP( 1, l)*l SORT! (X-XSEPI /Ix-XATT+RD.T) ) T 1.) 

RETURN 

END ‘ 
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* 

FUNCTION GAHSUMt Xt GAHA W» XIW #N WAKE ) 1143 

DIMENSION GAMAWI IOOO), XIWl 1000) .1144 

DATA PI / 3.141W 1145 

SUM=0. 1146 

NWM*NWAKE-1 . 1147 

DO 5 J =>2i NMM .. 1148 

JP=J«-l 1149 

DL= X IW I J 1 -X * 1150 

5 SUM*SUM*(GAMAWtJWGAMAinxi-GAHAHrjPl»*DL/tXn<IJPI-XIW<J),)l*ALOGI 1151 

1 (XIWI JPJ-XI/DU 1152 

ARG*l.-X .1153 

X M= ARG . 1154 A 

IF( ARG.LE.O.I ARC* 1 • . ' < 1155 

XP=l.+< 1 1 ' ... 1156 

AP*XP , . , 1157 

IFIAP.IE.O.) AP*1. . 1158 

D2=XIW< 2)-X ,1159 

GAMSUM*SUM+XM*(GAMAW{ 1 1-GAMAW I 2 ) » *ALOG I D2/ARGI / I.XI WI2 1^1. 1 1160 

1 ♦. 25*GAM AWIl )*t"4.*7M.OGtD27*XP* 1 1.-3. * X) *ALOG I API*6. *X ,1161 

1 -XM*(5.+3.*X)*AL0G<ARG)) 1162 

6AMSUM».3*6 A HS U H / P T 1163 

RETURN 1164 

END 1165 
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* 

SUBROUTINE W ASH2 < VZ , XGAM, NGAM ,C R , THl C K ,NF , ACAP. G AMI , I SE P , DE IP I , , 1166 

1 X I 2 ) 1167 

, DIMENS ION VZ ( 30), XGAM! 30),. THICM24) ,ACAP(30, 3) ,SN(30) ,CS(30) 1168 

DATA. NT Of 73/ - - , ■ . ' " 1169 

NTOP=NGAM , ... ■ ... , ...» . H70 

I F C ISEP .EQ.l) NTOP=NTOT ■ ■ .„ . H71 

t VZ( 1) = DEIP I*( .806853+AL0GIXI2-1.) ) „ ■- - 1172 

NGP1 = NGAM *1 , . .. • ... - . .. ,. . 1173 

S IGN= 1. ' . v 1174 

SUM=0 • ' • - 1 175 

. 00 5 N= 1, NF . . ' 1176 

SIGN=~SIGN ‘’'1177 

5 SUM=SUM+S IGN*N*THICK(N ) ; H 7 8 

F ACT=2. *< CR+2.*SUM ) j H 79 

SUM=0. • \ ,i 1160 

, . S IGN= 1. 1 H81 ' 

00 6 N* 1, NTOP | • . . \ H82 

S IGN=-S IGN 1183 

6 SUM=SUM+S IGN*N*ACAP (N*l, 1 ) . ;!• H84- 

VZ ( NGP l )=2.*CR*( SUM+ACAR (1,1)) *OELPI ♦ ( (XI 2* l.) *ALOG< ( X 12 +1 . )*. 5 )/ 1185 

1IXI2-1.)-1.5) ' . 1186 

VZ(NGP1) = VZ(NGP1H-FACT*ACAP(1,1) 1187 

10 NSUM=NGAM , - ■ i. 1188 

I F(NF .GT .NGAM ) NSUM=NF ... 1189 

DO 50 M = 2 * NGAM . ' 7 1190 

COST=XGAM(M) 1191 

T HET A® AR CO SI CO ST I i . ' . . r 1192 

DO 15 N = l , NSUM U 93 

ANGLE = N*T HET A 1194 

SN(N) = S IN (ANGLE) •. - 1195 

15 CS(N) = COS(ANGLE ) • 1196 

. S I NT= SN ( 1 ) 1197 

. C OT T = CD ST / 5 IN T 1198 

ST=0. . i ■ 1199 

S T P R= 0 . 1200 

SA=0. 1201 

SAPR=0. 1202 - 

DO 16 N= l , NF 1203 

ST=ST ♦THICK! N )*SN! N) 1204 

16 STPR=STPR*THICK(N)«‘(COTT*SN(N)*N*CS!N)) 1205 

DO 17 N=1,NT0P 1206 

NP1=N*1 1207 

SA=SA+ACAP(NPl, 1>*SN(N) 1208 

17 S APR=S APR +N*ACAP (NP1,1)*CS(N) 1209 

H=. 25*GAM 1*1 r. ♦COST l*r3.»CCTST-r.1- 1210 

HPR=-GAM1*SINT*( l.*1.5*C0ST) 1211 

F ACT 1=CR+STPR 1212 

FACT 2 = CR*( 1 CO ST ) •*■ SIN T* ST 1213 

T ERM1=ACAP( 1, 1)*( l.-COSTJ+SIMTMH+SA) 1214 

T ERM2= ACAP ( 1, 1 ) +COTT* ( H*-SA) +HPR* SAPR 1215 

40 XP=l.+COST 1216 

X M= 1 .-COST 1217 

50 VZ(M)=FACTl*TERMl+FACT2*TERM2+Otci'»*i i.3*^UST^.25*C1.-3.*C0ST)4XP* 1218 

1AL0GI XP/XM )♦( XI 2-C0ST )*ALOG! ( XI 2-COST )/ XM) /( XI 2-1 .) ) 1219 


RETURN 

END 
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* 

SUBROur INE UNPOPI NGP l f AR, AL AM ,AFACT,RMAT,CMAT,XGAM,AS ,MX,NZ,XC, 1222 

1 U£,NDIMC » AFAC2* AR2 ) 1223 

■ f COMMON /LOADS/ CR , CRHAT, THICK r24r*THATT24> ,NF ,ACAP 130 ,3 ) »ACAP2 (30 , 1224 

13 l , NG AM, GAM AW I 1000 I , GAMA W 21 1000 ) , XI Wt 1000) «NWAKE » XSEP » XATT , BCAP (10 1225 

10, 3), NS IG,R0B8,AA,BB,AKK, SS ,S SLAM ,2S ,UT , UN, UI NF ,XSIG(100) ,XS IGA (10 1226 

10 ),XS IGB( 100),NSIGA,NSIGB,DXI 1227 

DIMENSION VZ( 30J,AR2I30) 1228 

DIMENSION AR I 30 ), ALAMI 30) ,XGAM( 30), AS (3 0,3 01 , XC (300 I «UE(300,3 I 1229 

DOUBLE PRECISION RM AT ( 1 30 ) ,CHAT I NOIMC.NDI MC) 1230 

DO 5 M* l, N GP 1 1231 

SUB=ARIM)-ALAM(M)*AFACT/3. 1232 

RMAT ( M ) =SUB 1233 

CMATIM, 1)*1. , - 1234 

CMAT ( M , 2 )*XGAM ( M J 1235 

DO 5 N*2, NGAM ' 1236 

5 CMATIM, N+1)*ASIM,N) 1237 

CALL ALSOLI NGP It CMAT, RHaT, NDIMCI 1238 

DO 10 N® 1 » NGP l 1239 

10 ACAPIN, U-RMATINI 1240 

CALL WASH2(VZ,XGAM,NGAM,CR, THICK, NF,ACAP,0-,0,0., 2.) 1241 

DO 16 N * 1, NGP1 1242 

CMATIM, 11*1. 1243 

CMATIM, 2) * XGAMIM) 1244 

RMAT I M I * VZIMI— ALAMIM l*AFAC 2/ 3.+AR21M) 1245 

DO 16 N * 3, NGP1 1246 

16 CMATIM, N ) - ASIM, N-1I 1247 

CALL ALSOLlNGPirTHTWl RTTATi NOTHCI 1248 

DO 17 M* l , NGP 1 1249 

17 ACAP2I M, D-RMATCH) 1250 

GAM All I 1)*0. 1251 

G AMAW 21 1 ) *0. 1252 

DO 15 M* 1 » MX 1253 

SIGN»1. 1254 

IFIM-NZI 12,14,14 1255 

12 SIGN-SIGN 12 56 

14 CALL QECALtUE(M,l>,XCtM),SIGN?0) 1257 

15 CONT IN JE 1258 

RETURN 1259 

END 1260 


i 
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SUBROUTINE PR IMESI CP. TP, THETA ,CR ,THIC K.C AMB R» NF > 

DIMENSION THICK! 241. CAMBR«24» 

S ST*0 . 

SCT=0. 

SSR=0. 

S CR=0. 

ANGLE=0. 

C0UNT=0. 

DO 5 N= I, NF 

ANGLE=ANGLE+THETA I 

COUNT = COUNT *•! •* 

CN= COUNT *COS( ANGLE) 

SN=SIN( ANGLE) 

SST=SST ♦T HICK ! N )*SN 
SCT=SCT+THICK! N)*CN 
SSR=SSR«-CAMBR(N)*SN 
5 SCR=SCR*CAMBR(N)*CN 
C0ST=C3S( THETA) 

S INT=SIN( THETA) 

T P=CR*( C0ST-C0ST**2«-SINT**2HS1NT*(2.*C0ST*SST*SINT*SCT) 

C P= COS T * S S R *S INT*SCR 

RETURN 

END 
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* 

SUBROUTINE ATTPR I PR EC , ASZ , A S, AR ,C MAT, RMA.T , NG PI , AR2 , NDI MC , XGAM ,UTU , 1284 

1 PRECS ) .!• . . . 1285 

COMMON /LOADS/ CR.CRHAT, THICK ( 24) .THAT < 24) ,NF .ACAPI30 ,3 ) , ACAP2 130 , 1286 

13 1, NGAM, GAM AM I 1000 ) * GAMA W2( 10 00) * XI Ml 1000) , NWAKE »XSEP»XATT iBCAPHO 1287 
.10,3), NS IG,RDBB,AA,BB,AKK, SS,S$LAM,2S,UT,UN,UINF .XSIGI100I .XSIGAI10 1288 

.io I.KSIGBl 100).NSIGA,NSlGS,0xl 1289 

DIMENSION XGAMI30), A SZ< 30) ,A SI 3 0 ,30) , AR { 30) . V2 130) , AR2 (30 ) 1290 

DOUBLE PRECISION RMATI130), C MA T{ NO I MC , NDI MC ) 1291 

OATA PI /3. 14159/ 1292 

S AvE=XS IGINS IG+l) - 1293 

- XATS=XSIGA(NSIGA*1) 1294 

IEISAVE.LT.XATS) GO TO 20 1295 

PREC=PRECS*( l.-SAVE)/! l.-XATS) 1296 

RETURN • - 1297 

20 DO 50 H=1,NGP1 . 1298 

CMATIM, 1)=ASZ(M) 1299 

RMAT(M)=AR(M) r 1300 

DO 25 N= 1 , NGAM ; '• . ■ T 1301 

25 CMATIM, N+1) 3 ASIM,N) i 1302 

50 CONTINUE - . . •• . , • 1303 

CALL ALSOLHslGPir CHTTY RMAT, NDTHCT 1304 

DO 75 M 3 1 , NGP 1 1305 

75 ACAPIM, 1) 3 RMATIM) 1306 

GAMAHI 1)»GAM1(ACAP,0XI,PI ) 1307 

DOEL = UTU*GAMAW I 1) 1308 

DEL Pl 3 DDEL /P I 1309 

XSIGINSIG+D-2. 1310 

CALL MASH2I VZ, XGAM, NGAM.CR, THICK, NF, AC AP.GAMAWll) ,0,DELPI ,XIHI2) I 1311 

DO 60 M * 1, NGP 1 1312 

CMATIM, 1) * ASZ I M ) 1313 

RMATIM) 3 VZ IM ) EAR 21 H ) 1314 

DO 60 N * 2, NGP 1 1315 

60 cmatim, Nr* Trenrrw-Ti 1316 

CALL ALSOLt NGP 1, CM AT, RMAT, NDIMCI 1317 

DO 65 M*l, NGP 1 1318 

65 ACAP2IM, 1 )*RMAT(M> 1319 

GAMAM2I1) * GAM1IACAP2, DXI, PIU-DDEL 1320 

CALL CPCI PREC, SAVE, l.,0) 1321 

XSIGINS IG+I J'^SAVE ‘ 1322 

RETURN 1323 

END 1324 
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SUBROUTINE T AYLOR ( X SE P ,UE , XC . X , NZ , MX) 1325 

DIMENSION UE( 300,3), XC ( 300). XI300) ,, . • .1326 

DO 3 M = NZ « MX ' ~ '/•'”••• •• v .-1327 

' IFIXSEP.LT .XC1 M ) ) GO TO 4 ' "l328 

3 CONTINUE ■ Hi 329 

4 MS = M- 1 -1330 

OS EP=U E ( M i 1) -1331 

X STOP = X ( M S ) 13 32 

X BS EP = X STOP*! XSEP-XC! MS ) )* ! XI MI-XSTOP ) / ! XC !M) -XC ( MS I ) - 1333 

QMAX=0. 1334 

DO 5 M= NZ i MS 1335 

IF(UE(M, l) .LT.QMAX) GO TO 5 "1336 

QMAX=UE(M, 1) ' 1337 

MQM AX = M 1338 

5 CONTINUE 1339 

S M AX = 0 • - , • • 1340 

00 6 M=MQMAX , M S • . i i * IS - 1341 

M P= M+- 1 1342 

SLOPE=(UE(M,l)-UE«MP,in/(X(MPI-X(M)I 1343 

IF( SLOPE. LT.SMAX I GO TO 6 . 1344 

S MAX* SLOPE 1345 

MS M AX =M i 1346 

6 CONTINJE . t * ■ 1347 

IFIMSMAX.LT. MS ) GO TO 61 f 1348 

WR ITE( 6, 85 I XC ( MS ) 1349 

85 FORMAT! //10X, 'POSITION OF MAX SLOPE DOWNSTREAM OF X =*,E13.5,* SO 1350 

1N0 SMOOTHING NEEDED*) ■ 1351 

RETURN 1352 

61 MSM=MS- 1 1353 

DO 7 M=MSMAX,MSM 1354 

SLOPS=SLOPF 1355 

MP=M+1 1 1356 

SLOPE* I UE(H» lT-UErKP»T) T/TX IMP) — XTMT) 1357 

TLAM=( JE(M. H-QSEP l/l XSTOP-X(M) ) 1358 

I F( TL AM .GT .SLOPE I GO TO 8 1359 

7 CONTINJE 1360 

WRITE(6»88) , 1361 

88 FORMAT! //10X,' SLOPE MONOTONIC - NO SMOOTHING NEEDED*) 1362 

RETURN . • 1363 

8 SLOPE=SLOPS 1364 

K=M-1 1365 

K P=M 1366 

XK P=X ( <P ) 1367 

X IS=XBSEP-XKP 1368 

QKP=UE1 KP-, 11 1369 

QDI FF= QKP-QS EP 1370 

8=2 . *QDIFF/X IS— SLOPE 1371 

K PP=K P* 1 1372 

DO 9 M=KPP,MS 1373 

X I=XIMI- XKP. 1374 

WRITEI6.99) XC! KP ) t XC ! MS I 1376 

99 FORMAT ( / / 10X»"*EtOW "SMOOTHED ' BETWEEN - “X"** »E13.5 » * AND*.E13.5/I 1377 

RETURN 1378 


END 
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* 

SUBROUTINE ELDER! BCAP » XSI G»NS IG , UI NF , ELD »Y» YMAX , AZ ) 1380 

OIMENSION 8CAP( 100, 3) ,XSIG! 100) 1381 

DATA PI/ 3. 1*15926/ 1382 

BCAP1NS IG + 1, U»0. ' 1383 

XS=XSIG!l> 1384 

XZ*XSIG!N$ IG*l J 1385 

IF1XZ-1 • 1 16,16,1 1386 

1 OEADL=Xl-XS < . 138T 

NSIGM* NSIG-1 1388 

DEADL2* DEADL / 2. 1389 

Y MAX* 1 . E- 10 1390 

Y = PI ♦ DEADL 2 *1 aCAP(Ul) *■ BCAP 12 til / 2 .1 1391 

DO 30 J* 2, 1001 1392 

T* PI/1000. * (J-l) 1393 

YX= 8eAP( 1* 1) ♦! Pt-T * SIN1T1- 1 1394 

1 ♦ .5 * BCAP (2,11 *1 PI -T *.5* SIN! T*2. ) I 1395 

SUM*0 . 1396 

DO 31 I- 2, NSIGM 1397 

31 SUM * SUM ♦ BCAP! 1*1, 1 ) ♦ t SIN! T)/!I-1) 1398 

1 - SIN! (1*1) ♦TI/II41) 1 1399 

YX* DEADL2-*TYX-.5*SU«7 1400 

IF! YX .GT. YMAX I YMAX*YX 1401 

30 CONTINUE 1402 

ELD*Y/YMAX . 1403 

IF1ABS1EL0I-UINF) 2O,20r,lZ 1404 

12 IF! ELD) 14, 16, 16 1405 

14 ELD“-UINF' 1406 

GO TO 20 1407 

16 ELD=U1NF ' / 1408 

20 CONTINUE 7 1409 

MRITE1 6, 801 EL CT 1410 

80 FORMAT! /10X, 7HELD0T -E13.5/I ./ 1411 

RETURN- •' ’ 1412 

END ' ; ' 1413 
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SUBROUTINE R EATTI UC, V, X, Y,M X,N Y ,R Y ,DR Y,UE ,X5 ,DE L5 ,MST ,REB» , 
DIMENSION UC( 100, 3), V( 100,21, Yl 1001 
DIMENSION XI 300T,UEC 300, 31 

0 I MENS ION 'TABU 24) # TAB 2( 24) # TAB 3 1 24) # TAB4 (24) • XI TAB (24) 

DATA TAB1 /24.9S, 23.29,21 . 04, 19. 33# 17.61 ,15.29 ,13.46,1 1.54, 1 0.36*9 

1.38,8.35,7.32,6.29, 5.31,4.4,3.57 ,2.22 ,1.26,. 66, .31 ,.14 ,.01 ,0. ,0./ 
DATA TAB2 /20.05, 18.85,17.25, 16.04,14.8,13.12,11.77,10.3,9.36,8.65 


1,7.95,7.2, 6.43, 5.66,4.9,4.18, 2. 89, 1. 86,1 .1 1 ,. 62 ,.32 ,.04,0. ,0*/ . 
DATA TAB3 716.65, 15.8, IT. 67,13.8,12.91, 11. 66,10.65 ,9.48 ,8.71 ,8 .1 1 , 

17.59.7.01. 6.41. 5.77. 5. 13.4.5.3.31 .2. 2 8. 1.48.. 9.. 51 .. 09.. 01.0./ 

DATA TAB4 / io . 12, 10 .05, 9. 93, 9. 78, 9. 58 ,9. 17 ,8. 72 ,8. 08 ,7.6 ,7.2 ,6.85 , 

16.53.6.18. 5. 79. 5. 36. 4. 91. 3. 98. 3. 05. 2. 21. 1.5.. 95 ..22 ..03.0.7 

DATA XITAB /.0001, .0002,. 0005, .001, .002, .005, .01, .02, .03, .04* .05,. 
10 6, . 07, • 08, . 09, . 1, . 12 , . 14 , . 16 , . 1 8 , . 2 , . 25 , . 3 , . 3 5 / , 

3 FORMAT ( ///40X.23HAT REATTACHMENT, BETA *E13.5) 

MOUT=6 •■'''■ 

’ y ' RTR=SQRT(REB) 

UC ( 1* 2 )=0 . 

UC! 1, 3)*0 . 

V ( 1, 1 ) *0. 

VII, 21*0. 

DO 5 M* 1, MX 


I F( X5-X (MI) 4,4,5 

4 MST=>M ♦2 
GO TO 6 

5 CONTINUE 

6 X A«X(MST-2J 
X B*X I MST-1 ) 

U A*UE{ MST-2, 1) 

U 8*UE! MST-1, 1 1 
Z A* AL0G1 U A*DEL 5*REB ) 

PGRAD*2.*( UA-UBl/l ( UA+UB )*( XB-XA) 1 

BETH2*! . 0974- SORT! DEL 5*P GR A 0 J )/ < . 0249+i 004565*ZA» 


I F( BET M2- 1. ) 8,7,7 

7 BET M2* 1 . 

GO TO 10 

8 I F( BET M2- .31 9,9,10 

9 • BETM2* • 3 

10 BETA*l./t BET*2*BFnr2) 



WRITE! NOUT, 31 BETA 

AGAM* .09?4*BETM2- . 0249/BE TA 

BGAM*. 304565/BETA 

AH*l.-( 5 »3+3v9*1TETN~2I*t ,0974^ . 0249*8 E TM2 ) 


BH*BETM2*( 5.3* , 3.9*8ETM2)* .004565 
GAMA* AGAMF BG AM *Z A 

DER IV=J A*REB*EXP( -ZA »*GAMA* GA MA* < l.+BE TA*( 1. *AH«-BH*ZA> > / 1 AH*8H*BH* 


1ZAI 

ZB*ZA+DERIV*!XB-XA1 
DEL 8= EXP (ZB) /TUB^REBT 
GAMB*AGAM-BGAM*ZB 


11 I FI DELL-Y ( NY-3 ) l 14,12,12 

12 RY*RY«-DRY 

CALL YSET(RY,Y(21»NV,Y> 


* 
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GO TO 11 

14 I FI BET A-4 • ) 102,101,101 

101 TERPB=l.-4./8ETA 
INDEX = 3 

GO TO 110 

102 IF( BETA-2.) 104, 103,103 

103 TERPB=.5*BETA-1- 
!N0EX=2 

GO TO 110 

104 TERPB*8ET«-t. 

I N0EX=1 

110 K=0 

T ERP 1* l ,-T ERPB 
50 K=K+1 

GO TO C 16, 17, 99), K 

16 G=GAH A 
DELT A=DEL 5 
l/EOGE*JA 
L= 3 

GO TO 18 

17 G=GAHB 
DELTA* DH.B 
U EDGE* JB 
1*2 

18 X IC0*G/(DELTA*RTR*BETM2) 
U COW* RT R« ( UEOGE*G T»» 2 
EFC0=G/BETM2 


NLAK*NY 
DO 75 M*2,NY 
XI-VCNI4XIC0 - 
IFIXI-. 35) 20, 19,19 
19 UC(Na)*t#EOG€ 


GO TO 75 
20 CALL 

INDPl* INDEX -fl 
CALL 

FP*TERP l*FPl*TERPB*FP 2 


t,F PIT 
P2> 


UCIN,rT«UED G E»t l.-EF CP »FP7 
IFCN-N LAM I 21, 7 5, 75 

I FI ALTER-UCIN.L)) 33,33,32 

32 JC(N,T7*AtTB. 

GO TO 75 

33 NLAH*K 
75 C0NT1NJE 

GfT TD-50 — 

99 DO 60 K*2, 3 


S AVE2=0 - 
DO 60 N*3,NY 
SAVEl*iCtN- UK ) - 
UC(N-1,K)»(SAVE2*SAVEUUC(N,IU ) /3 


60 


DUDX*0. 

C0D*.5/TX«=**1 
00 65 N*2,NY 
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1500 
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OUOXP=COD*'{UC( N# 2 l-UCl N* 3 ) ) 

V(N»LI = V( N— 1 « 1 1-< V( NI-YCN-1 )) *( DUDXP+DUDXJ 

V { N» 2 ) = V ( N* 11 

DUDX= OJ DXP 

RETURN 

END 


1523 

1524 

1525 

1526 
152? 
1528 



* 

SUBROUT INE ELP TT1 ALPH1 , ALPH2, EH l ,TURF ,THET2 ,Ul NF ,0X1 ,CHPA,CHPAS» 1529 

S AV ET = AL P H 1 1530 

STEP*T3RF*DXI 1531 

SINS*SIN<STEPI 1532 

COSS=C3S« STEP ) 1533 

CONST “2 .*EM 1*1 UINF/TORF 1**2 1534 

ALPHl = T HETZ ♦( ALPH1-THET2 1*C0S S+ALPH2*3I NS/T0RF*C0NST*(2.*CMP4-CHPA 1535 

IS )*( 1.-C0SS)+C0NSTP(CMPAS-CMPA1*1$1NS-STEP*C0SSI/«T0RF*DXI) 1536 

ALPH2 = ALPH2*C05S=TtJRF*SlNS*rSAVET-THETZr*C0NST* CCHPA-CMPAS1*! l.-CC 1537 

1SS)/DXI*C0NST*CNPA*T0RF*SINS 1538 

RETURN 1539 

END 1540 


83 



* 

SUBROUTINE VWASHI BAR G* H* S ,N VOR , XI ,UI NF * VZI P»XGAK,NG PI t OXI ) 1541 

DIMENSION VZ I P C 30)»XGAMt30) 1542 

DO 10 N = 1 # NGP 1 1543 

DIFF=XGAM(N)-X1 - 1544 

SUM=0 . 1545 

•* DO 5 K= 1» NVOR 1546 

SUM=SUM+D IFF/< OIFF*DIFF*H ) 1547 

5 DIFF*DIFF-S 1548 

10 VZ IP< N) =VZ IPIN M-SUM*BARG 1549 

RETURN 1550 

, END 1551 
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• ; . • . * 

SUBROUTINE W AS H< XGAM, NGA M , T I ME , A LPH1 , ALPH2 .HEAVE , ARCT ,FRECF ,PHIH ,U 1552 
11NF,CA'1BR,NF,VZIP,M0TR,IN0V) 1553 

, DIMENSION XGAMI 30J.VZ IPT JDT.CAMBR (241 1554 

NGP1 = NGAM+1 1555 

ANGLE = FREQF*Tt>f£ 1556 

GO TO (108.120). INOV 1557 

108 j. GO TO ( 110, 120I.M0TR 1558 

110 CONST =-ALPH2*C0S( ANGLE )*UINF +HEA VE*CUS( ANGLE^PHI H ) «-ALPHl*U INF . 1559 

FACT =-ALPX2*FRE0T r *SIN(ANGLFr*UINF 1560 

- GO TO 130 1561 

120 CONST=J INF *AL'PH1+ HEAVE 1562 

FACT=-J 1NF4ALPH2 1563 

130 DO 10 1 , NGP 1 1564 

X=X GAM CM) 1565 

THETA = ARCTIX) 1566 

SUM=0. 1567 

COUNT = 0 • 1568 

DO 20 N* 1 , NF 1569 

COUNT = COUNT ♦ li 1570 

20 SUM*SU1+C0UNT *CAMBR( N )*COS( COUNT* THE TA) 1571 

IF(M-l) ^;4f2 1572 

2 IF(NGPl-M) 3,4,3 1573 

4 SUM = SUM V SUM 1574 

GO TO 50 1575 

3 COUNT * 0. 1576 

COTT = X/S IN(THETA ) 1577 

, DO 30 N-WNF-" 1578 

COUNT - COUNT+THETA 1579 

30 SUM=SUH*COTT»CAMBRTWT*S IN t COUNT » 1580 

50 VZIP(M) = UINF*SUM+C0NST+FACT*(AR0T-X) 1581 

10 CONTINUE 1582 

RETURN 1583 

END 1584 
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SUBROUT INE M IX ERC FPRE S ,PR EC , UINF , UDOT.THI CK ,NF , XBSI G , NS IG , I NOT, DEL 
lltTHET 1,REB,USEP,X4,CP4, ROBB) 

DIMENSION FPRESMOO ) » THICK! 24 ) , XB SI G 1 1 00) 

FCAPI X )=-19.556*'X + 107.535«‘X*X-336.33*X**3*508.l*X**4-295.96*X**5 
0 IUX J = -.46532*X«-.68425*X*X-.45293*X**3*.6592*X**4 
U I2IX I=-.045929*X-1.91615*X*X*2.91843*X**3-5.42125*X**4 
CR=SQRTIR0B8)/4. 

DIST*=.5*I XBS IGI 2I-XBSIGI 1> ) 

XSEP=XBS IG( D-DIST , 

X ATT=XBSIG(:NSIG)*PIST 

IF INDT IS NONZERO, THE BOUNDARY LAYER IS TURBULENT 
AT S EPARAT ION. \ 

CALL H4X4 ( IN DT ,XSEP,DEL1, THET 1 , X ATT, RE B , USEP , X3 »H3 » X4 »H4) 
IFIXSEP-I.) 24,25,25 
2 ; 5 CP4=0. i : . 

GO TO 27 - 

24 URAT=EX PI -.08712-UI II H4I- .24723* I. 32 55+ UI 2 IH4 ) ) ] ,, 

CP4= l.-( l.-PR£C)/URAT**2 
DEADL=XATT-XSEP 
IF( DEADL-2.I. 5,6,6 \ 

5 G=( .5*DEADL 1**2 ; 

GO TO 7 « 

6 G= l . . : 

7 CP4=PREC + (CP4-PRECI*I l.-G*XSEP) 

27 CONTINJE 

CP4L IN = 2. ♦ I 1.-SQRTII.-CP4I » .. 

PRL IN = 2 .*( 1.- SORT l l.-PREC ) I 
COE F= ( PRL IN-CP4L IN 1/ IXATT- X4I 
CZ»2.*J00T/UINF ' , \ . 

C2 = -2.*U INF -/ . r . ' 

DO 20 M* 1, NS IG 
SUM=0. 

X=XBSIG(M) 

IFIX-1.) 2,2,3 

2 THETA = ARCTIXJ : 5 

COST=COS( THETA) - 

S INT=S INI THET A I ~ , ,, 

COUNT *0 . Y 

DO 10 N=I»NF . t 

C0UNT=C0UNT*1. 

ANGLE=COUNT*THETA 

10 - SUM=SUM-*-T HICK! N )♦( COUNT* S IN T* SIN I ANGLE )-2.*C0ST *COS I ANGLE I I 
SUM=-4.*U INF»T SUfTfCR^l.^Ti^OTST) ; 

GO TO 35 

3 SQSS=SORT I X*X- 1. ) 

XRAO=i./I X+SQSS) 

COUNT = 1 . ; >' ' ' 

fact-xrad / ' i..‘i _ 

DO 30 N=2* NF 
count* count - n. 

F ACT= FACT *XRAD 

30 SUHaSlM 4-THICKIN )*FACT*(2.*X— C OUNT* SQS S) 
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SUM=SUH*2.*2.*CR*<X*X- l.-SQSS*( l. + X+XJ / 1 l.+XJ) «-3.*THI CK(1 ) *XR A0**2 1639 

$UM=5 JH*U INF 1640 

35 CP=CP4L IN 1641 

I Ft X-X 4 J 55.50.'5*r~ 1642 

50 CP«CP+< X-X4I*C0EF 1643 

55 CONTINUE 1644 

FPRESI m»-UINF*CP*SUH 1645 

20 CONTINJE 1646 

RETURN 1647 

END 1648 
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* 

SUBROUT IN E SECKXU, YU, XL . YL ,NOFF , NF ,RDBC , ST, SCI 1649 

C PROGRAM TO COMPUTE COEFFICIENTS TN AND CN OF THE FOURIER SERIES 1650 

C REPRESENTATION OF SECTION THICKNESS AND CAMBER DISTRIBUTIONS. 1651 

DIMENSION XUI 30 ) , YU l 30 ) , XU 30 ) , YL I 30) .YUC130) .YLCI30) ,ST124) ,SC124 1652 

1 ) ,DUMl 50 ) ,TBAR( 50),CBARl 50) 1653 

11 FORMATl 1HI,4X*6HRDBC =E 13.5/4X, 7HRCDBC =E13.5) 1654 

12 FORMATl ////47X, 26HINPUT AND COMPUTED OFFSETS/) 1655 

13 FORMATl 19X, 4HX 1/C, 12X, 4HYU/C . 11X.5HYUC/C ,2 OX ,4H XI /C , 1 2 X ,4H YL/ C ,1 IX 1656 

1,5HYL C/C/) 1657 

14 FORMAT I 8X , 3F 16. 5, 8X, 3F 16. 5) 1658 

N A=6 1659 

RMA=6. 1660 

RNF=FL3AT INF) 1661 

MOUT =6 1662 

PI=3. 14159 1663 

D£LT=P 1/ ( 2.*RNF) 1664 

NT C=2*NF- 1 1665 

N INT=NTC+2 1666 

NS I MP=NTC +1 1667 

T HET A=0 . 1668 

DO 89 K=1,NTC 1669 

THETA=THETA+DElT 1670 

X 1=.5*( 1. +0051 THETA) I 1671 

00 90 LAM=2,N0FF 1672 

1 FC Xl-XUl LAM ) ) 110,90,90 , 1673 

110 YUINT=YU( LAM— 1 ) ♦( X 1— Xul L AM— 1) ) * I YUI LAM) — YUILAM— 1 ) ) / IXU(LAM)— XUILAM 1674 

1-1)) 1675 

GO TO 111 1676 

90 CONTINJE 1677 

111 DO 80 L AM=2, NOFF 1678 

IFIXl-XLlLAM)) 210,80,80 1679 

210 Y L INT=YL ( LAM— 1 )♦( Xl-XL(LAM-l) ) * I YL { LAM) -YL <L AM-1 ) ) / ( XL C LAM) -X L l LAM 1680 

1-1)) 1681 

GO TO 112 1682 

80 CONTINJE 1683 

112 TBAR(K*1)=.5*(YUINT-YLINT> 1684 

89 CBARI K «■ 1 ) = .5* I YU INT+YL INT ) 1685 

T BARI l ) =0 . 1686 

CBARI i)=0. 1687 

TBARININT )=0. 1688 

CBARI NINT )=0. 1689 

S AVT 2=0 . 1690 

S AVC2 = D • 1691 

DO 39 1 = 2, NS IMP 1692 

SAVTl = SAVrr 1693 

S AVCl = S AV C2 1694 

S AVT2 = T BARI I ) 1695 

S AVC2=C8ARl l ) 1696 

T BARI I ) = I SAVT 1+SAVT 2+TBAR 1 1 ♦T) 1/3. 1697 

39 CBARIJ) = ISAVCl+SAVC2«CBARlI*i))/3. 1698 

T T A=TBARl NA1 1699 

T T B=T BARI NA-*- 1 ) 1700 

- TT C=T BARI N A+21 1 7 °1 

T AA=DELT *1 RNA- 1. ) 1702 


456 


T BB=T AA ♦CELT , 

T CC=TBB *D ELT , 

X A=.5*C0SC TAA ) 

xB*.5*cos(T8irr , *• ' r. 

XC=.5*C0SCTCC> 

SLOPE =( (TT C-TTBT*( XB-XA) / ( XC- XB ) ♦CTTB-TTA) * l XC- XB) / C XB— XA ) ) / I XC-X A 
1 ) 

THETA=0. 

COSB=COS( TBB ) 

DO 456 I*2,NA“ 

T HET A=r HET A+0EL7 
COST* CO SC THETA) ~ 

T BARC I ) = ( SORT C l.-COST ) /1 1 .-CO SB ) **l . 5)* f TTB* (l.*COST-2. •COSBi/lli- 
icosBj+.5*SLOpE*tcosr-ct)SB r: 

T BARC I )*T8AR(I>*< l.-COST) 

N L E*2*M F* t-NA 
TT A=T BARC NLE+1 ) 

TT B*T BARC NLE ) 

TTC*TBAR< NLE-1) 

T AA*P I-DELT*C RNA— 1. ) 

T BB*T AA-DELT , „ , 

TCC=TBB-0€LT 
X A*.5*C l.+COSC TAA)) 

XB=.5*{ l.«tOSTTBBT) 

XC=.5*C l.+COSC TCC)) ... 

S L Q PE* C < TTC-TTBt»nC8- XAI 7 CXC- XB^M^TTB - TT A ) * | XC - XB ) / C X B- X A ) ) / 1 XC-X A 
1 ) 
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CAPA=SCOPF*SCJRTI 2.*XB/R0BCT 
C APB*TT B/SORT ( 2«*XB*RDBC ) 

AS= C 2 . 5*CAPB- .5* CAPA- 2.1 / XB 
BS=C .5*CAPA-1.5*CAPB*1.»/XB**2 
THETASPI — 

COEF=SQRT ( 2.PRD8C ) 

DO 457~ T=27Nt 
lND=2*MF+2-I 
T HET A=THET A-OELT 
X= . 5*C l.*COS< THETA) ) 

457 TBARI 1*0 )=COtF*SQimx^*t t^*AS***8 $**♦ X) 
COEF*COEF/4. 

T HETA-CK 

DO 458 I=2,N$I«P 
THETA*THETATtJECT- 
^ S INT=SINC THETA ) ^ 

T HET A*0 . 

DO 459 1^ 2, MS IMP 
THETA=THETA*OELT 
459 CBARC I ) *CWirtTTrStHtTH€TA )~ 

RKK*0. 

DO 59 X * 1 » NF 
RKK*RK<>1. 

T HET A=0. — 

DO 777 !»1, NINT 

777 ThETA=THETA*OELT 
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CALL SIMP(NSIMP,DELT,DUM, VARY) 1758 

S T ( K ) = 2 » *V ARY /PI 1759 

THETA = 0. 1760 

DO 888 1= 1»N INT 1761 

DUM( I ) = CBAR( I )*SIN( THETA*RKK) 1762 

888 T MET A=f HET A+-OEL T 1763 

CALL SIMPtNSIMP, DEL T, OUM , VAR Y ) 1764 

59 SC(K)=2.*V ARY/PI 1765 

DO 969 I = 1» NOFF 1766 

X=XU(I» 1767 

CALL EV AL(NF,X,SC,ST,CB,TB,COEF ) 1768 

969 JtUC< I ) = CB+TB 1769 

DO 869 { = 1, NO FF 1770 

X=XL( I ) 1771 

CALL EYAL(NF,X f SCtST.CB.TB.COEF) 1772 

869 YLC( 1 )=CB-TB 1773 

RS AV = R DBC 1774 

WRITE( NOUT, 11) RSAV.RDBC 1775 

HRITEINOUT, 12) 1776 

WRITE< MOUT, 13) 1777 

HR ITE( NOUT , 14) ( XU( I ) * YU( I ) t YUC 1 1 ) r XL (I ) ,YL i I ) « YLC ( I ) ,1 = 1 , NOFF) 1778 

RETURN 1779 

END 1780 
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* 

SUBROUTINE EVALINNF* XX»SSC»SST ,CCB ,TTB»COEF) 1781 

0 I MENS ION SSCI 50), SST( 50) 1782 

C0$T=2.*XX-1. 1783 

C0STS=C0ST**2 1784 

l F( COSTS- l.E-8 ) 303, 304, 304 1785 

304 TANT = SQRT( l. /COSTS- 1. J 1786 

THE = AT ANl TANT) 1787 

GO TO 305 1788 

303 THE=1.5708 1789 

305 I FI COST ) 403,404,404 1790 

403 THE=3. 14159-THE 1791 

404 ARG=0. 1792 

SU*l=0. 1793 

SUM2=0. ' 1794 

00 551 N= 1,'NNF 1795 

ARG=ARG-*THE 1796 

SUM1«SUH1*SSCIN)*SINIARG) 1797 

551 SUM2=SUM24SST(N)*SINIARG) 1798 

S INT=SIN< THE ) 1799 

CCB*SUHi*SINT 1800 

T T B= S I NT* t CO* WT4“C0 S T ) ♦SIWT»S t ) H 2 » 1801 

RETURN 1802 

END 1803 
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* 

SUBROUTINE AL SOL I NT • C, R, NDIMC) ' ' , ’ 180* 

DOUBLE PRECISION C ( NO IMC ,NDI MC ) , RU30I ' . ,1805 

DOUBLE PRECISION CMAX, SAVE, SUIT ’ 1 i . , ' 1806 

NT 1 « NT-1 1807 

DO 99 J»1,NT1 1808 

CMAX * C(NT» J 1 1809 

L =NT 1810 

DO 10 l 3 J,NTl 1811 

IF (DABS( CHAX)-DABS( C( I,J J ) ) 5,10,10 1812 

5 CMAX 3 C< I, J ) 1813 

L 3 1 1814 

10 CONTINJE 1815 

DO 15 J J*J , NT 1816 

SAVE * CIL.JJI 1817 

C(L.JJ) 3 CfJ.JJJ 1818 

15 C(J,JJ) 3 SAVE/CMAX 1819 

SAVE 3 R ( L ) 1820 

R ( L ) 3 R(J] 1821 

R( J» 3 SAVE/CMAX 1822 

JP1 3 J+i 1823 

00 25 WPliNT 1824 

DO 20 J J*JP1,NT 1825 

20 C(I,JJ) 3 C(IfJJ) - Cl 1,J)*C( J,JJ) 1826 

25 R ( I ) 3 RID - RI J )*C( I , J ! 1827 

99 CONTINUE 1828 

R(NT) 3 R( NT )/CINT,NT ) 1829 

DO 150 K«I,NTr ~ 1830 

I 3 NT-K 1831 

IP1 3 1*1 1832 

SUM 3 D. 1833 

DO 12 5 J=rPl,NT 18 34 

125 SUM 3 SUM * R(J)*C(I,J) 1835 

150 R( I) 3 RCTr=“SW“ 1836 

RETURN 1837 

END 1838 
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* 

FUNCTION GAM1(ACAP,DXI,PI ) 1839 

DIMENSION ACAPI30,3i 1840 

GAM1»PI*( -I.5*TtC7tP(Ttn-.75*ACAPl2,ri^2.*ACAPIl f 2)*ACAP<2 ,2)-.5*AC 1841 

IAPI I, 3»-.25*ACAP( 2, 31 l/DXI 1842 

RETURN 1843 

END 1844 


i 
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• * 

SUBROUTINE EGAM II NU, NG , A , B, XSEP . XA TT, GAMMA , Y,GI I . 1845 

DIMENSION A( 30,31 18<,6 

S INT=SQRT ( 1 Y*Y ) 1847 

T HET A= ARCTI Y I 1848 

SUM=0. 1849 

• COUNT = 1 . 1850 

DO 6 N=2*NG ’ 1851 

C0UNT=C0UNT*1. 1852 

6 SUM=SUl4-AIN*l.NU)*(SIN<<C0UNT+l.)*THETA) /TC0UNT*1 . » -SI Nl ICOUKT-1 * ) 1853 

1*THETA) /( COUNT- I. )) 1854 

GI = (3.14l59-THETA«-SINT!*(Al 1 ,NU )♦. 5*A ( 2 , NU 1 ) ♦. 5 *SUM-.25*GAMMA* I l .♦ 1855 

IY ) *S I N T *S INT 1856 

I F ( Y-X ATT ) 8,8,7 1857 

7 DIFF=1.-XATT 1858 

IF( 0IFF-1.E-6) 8,8,9 1859 

9 „ GI = GH-2.*B40IFF**(- 1.5I*SQRT( ( XATT- XSEP1 *( U- Yl * ( Y-XATT >1 1860 

8 , CONTINJE 1861 

> RETURN 1862 

END 1863 
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* 

SUBROUT INE BUBB( DEL 1, THE T 1 ,RE B , XC 1 ,Ul , XC 5 ,DC P ,DE L5 , X , XC , MX ,NZ,X5,U 18 64 
15»UE,ALTC»RENEL*USTOP) 1865 

DIMENSION Xt30(jr.XtT3O0),0E!300,3) 1866 

FCAPIX )=- 19.556*X*107. 535*X*X-336. 33* X** 3*508 .1 *X**4-295. 96*X**5 1867 

U I1IX )=-.46532*X*.68425*X*X-.45293*X**3*.6592*X**4 1868 

U I 2<X ) =-• 045929*X- 1 .9 161 5*X*X*2 • 9 1 843* X* *3-5. 42125*X**4 1869 

FOELT! X ) = EXP< 2.5773-. 342 52*X- ,4379*X* X-. 07651 1* X**3-. 0039707 *X**4 ) 1870 

FAICHI X )* EXP( -3 . 74 81*. 038772* X*. 4 196 7*X* X*. 071046*X**3*.0032162*X* 18 71 

1*4) 1872 

OEL I(X )®-.045929*AL0GI X )-3. 9242*X*. 54535*X* X-l. 39147*X**3-10. 8425* 18 73 

IX **4 18 74 

25 FORMAT! IHl, 44X, 31HANALYSIS OF LEAO ING-EOGE BUBBLE ////34X« IHX 1 19X » 1 1875 

IHU « 19X , 1HH» 18X, 4HD I'SP / ) 18 76 

30 FORMAT! 20X.4E20. 5) 1877 

M0UT*6 1878 

Hi® .25 1879 

H5= .429 I860 

DO 5 M*NZ • MX 1881 

IFIXCI-XCIMH «-■*. 5 1882 

4 Ml=M 1883 

GO TO 6 - 1884 

5 CONTINUE 1885 

6 Xl*XlHl-l)*!xrHl)-XTHr-ri )*IXCl-XCim-ll »/!XC!Ml)-XC!Ml-l) ) 1886 

X4®Xl*RENEL/!Ul*REB) 1887 

ARG*AL3Gt LX*-XT)VtREB*DEL t*OELl*Ul)) 1888 

H4= .25*FAICHI ARC) 1889 

D EL 4* . 5 8*F DEL T I ARG)*DEL 1 1890 

X 5®X4* 10 • 5*0EL 4* ! 1.-! H4/. 4291**21 1891 

I F(U1-UST0P» 41*47.40 1892 

40 ALTL*ALTC*DELl 1893 

I F! X5-X 1 .LTJALTL'I" X5*X 1-fAtTt 1894 

41 URAT*EXP! - .087 12- U 1 1( H4 )- .247 23* ( • 32 55*UI2 !H4 ) ) ) 1895 

OCP=»UT*Ul*ri .-UR AT*»2T 1896 

DRAT® EXP ( -2.24374-FCAP1 H4 )♦• ,24723* (2 . 0214+DE LI IH4) I ) 1897 

DEL5«0RAT*0EL^' 1898 

DO 7 M»NZ »MX 1899 

I F! X5-XTM17 10,76,7 1900 

16 M5-M , 1901 

GO Tfr O 1902 

7 CONTINUE 1903 

8 FACT«!X5^7tTr5^trrAtXt1fST-XrM5^in 1904 

FACTl® 1 .- FACT 1905 

XC5»XC! M5-i ) »FA C T L l *X € !M5 ) » FA E7 1906 

U5*UE! 4 5-1, 1 )*FACT 1*UE !M5»1)*FACT 1907 

WRtTEt M OUT, 25) 1908 

WRITE! 40UT ,30) XI, Ul, HI, DEL 1 1909 

WRITE1M0UT i~307 X4,U1, H4i P6t'4 19 10 

WR ITE! 40UT » 30 ) X5,U5,H5,DEL5 1911 

RETURN 1912 

END , 1913 
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SUBROUTINE YSETl R . A ,NY ,Y ) . 

DIMENSION Y ( IOOI 
RP1=1.*R 

Y ( 1)*0. 

Y ( 2 )* A 

DO 10 N“3»NY 

YlN>«RPl*Y{N-i|-R*Y(N-2) 

RETURN 

END 


1914 

1915 
19 16 
19 17 

1918 

1919 

1920 

1921 

1922 



o o o O 


SUBROUT INE H4X4I INDT, X 1, DEL 1, THE T1 , X5 , REB ,U1 ,X3 ,H3 ,X4,H4> 

CURL FI H 1=26. 703/H+305.03*AL0GIH 1- 2111. 3*H*3327. 8*H*H-2403 .9*H**3 
FDELT I X )=EXPl 2.5773-.34Z52*X- .437 9* X* X-. 07651 1*X**3-.0039707*X**4> 
FA ICH( X J = EXP(-3. 7481*. 038772* X+ • 41967* X* X+. 07 1046*X**3*.00 321 62 *X* 
1*41 

10 FORMAT l //20X, 54HA SOLUTION FOR X4 COULD NOT BE OBTAINED IN 1000 TP 
HALS) 

M0UT=6 ’ 

;iF [NOT IS NONZERO, THE BOUNDARY LAYER IS TURBULENT 
AT SEPARATION. 


2 


5 


15 


20 
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41 

42 

43 


51 

52 
61 


50 


I FI INDT ) 2,5,2 
H3=THET 1/ DEL l 
X 3=X I 
DEL3=DEL 1 


GO TO 20 

X3*X1*5.E4/IU1*REB) 

AR G* AL 3 G( t X 3-X time B*D€ L I*DE L l ») 
H3=THET 1*FA ICHI ARG 1 /DEL 1 

IFIX3-X5) 20,15,15 

H4= .429 

X4“X5 

GO TO 50 

CONT IN J E 

I G0=0 

DIST=X5-Xi 
UNDER=0 . 

H4=»H3*H3 
COEF I® DEL 3*WJ 
C0EF2=10.5*DEL3*H3 
SUB*X3- C 0 EF 1* CU RL F ( H3 I 
0VER=»H4 


H4* .5*1 H4+UNOER 7 

X 4° CURL FI H4)*C0EFl*SUB 

ALT ER=X 5-C0ef ! ?*tt-.^tH4/i 42^7*42 J /K4 

IG0=IG3*1 

i n X4- ALTETT1 4T, 50,~4Z” 

IF! IGQ-1000) 95,61,61 

IFl ABSIX4-AL TER T/DTST- .OOI) 50, 50, 43 

UNDER=H4 


H4* . 5*t 0Vm*H47- 
X4=CUR.LFI H4)*C0EFl*SUB 


I GO* I GO *1 


I F I X 4- AL TER ) 52,50,^1 
IF! IG0-1000) 43,61,61 

I F I A B S I X 4- At TER ITtTLST- .001 I 50, 5 0 , 9 5 

H4* .429 

X4*rX5— 

WRITEI40UT, 101 


CONTINUE 

RETURN 


1923 

1924 

1925 

1926 

1927 

1928 

1929 

1930 

1931 
19 32 

1933 

1934 

1935 

1936 

1937 

1938 
19 39 

1940 

1941 

1942 

1943 

1944 

1945 

1946 

1947 

1 948 

1949 

1950 

1951 

1952 

1953 

1954 

1955 

1956 

1957 

1958 

1959 

1960 

1961 

1962 

1963 

1964 

1965 

1966 

1967 

1968 

1969 

1970 

1971 

1972 

1973 

1974 

1975 

1976 



, SUBROUTINE S ET SX( NSP 1 * XSEP * XA TT • XSIG (ANGLE ) 1978 

DIMENSION XSIGIIOO) 1979 

A* . 5* ( X S EP +XATT 1 I960 

B=.5*(XATT-XSEP> 1981 

ARG=0. 1982 

DO 5 N = 1 » N SP 1 1983 

XS IG(NI=A-8*COSC ARGI 1984 

ARG=ARG+ANGLE 1985 

RETURN 1986 

END 1987 
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* 


FUNCTION ARCTIX) 

1988 

PI= 3.14159 

1989 

I F ( ABS< X J-l.E-6) 1,2,2 

1990 

ARCT= . 5*P I 

1991 

GO TO 6 

1992 

IFIX+. 99999) 3,4,4 

1993 

ARCT=PI 

1994 

GO TO 6 

1995 

ARCT=AT AN( SQRT1 i.-X*XJ/X) 

1996 

IF(ARCT) 5,6,6 

1997 

APCT=ARCT *P I 

1998 

continue 

1999 

return 

2000 

END 

2001 
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SUBPOUT INE SCAL ! SBL * NSBL *F RZ » ARR * RDBB > 
OIMENSION SBL ( 300 ! 

D6LZ=FRZ *RDBB 
EN= ARR/ FRZ 
DO 5 N= 1, 300 
IHEN-N) 4,4,5 
N E=N 
GO TO 6 
CONTINUE 
NG= NS BL-NE 
EN= FLOAT! NG) 

NGM 1=NG- 1 

S BL ! 1 I = 0 . 

00 7 N=2, NE 

SBL!N) = SBL(N-1)+DELZ 

FP ACT =2 . 2/OELZ 

FRAC1=FRACT-1. 

R=FPACT**( 1. /FLOAT! NGMl) ) 

S AV E= R 

R = R-| R **NG-FRACT*R +FRAC 1 ) /! EN*R**NGH1-FRAC T) 

IF! ABS! S AVE-R I-l.E-6 1 9, 8 

RP1=R*1. 

DO 10 N=NE,NSBL 

SBL(N+1)=RP1*S8L(N) -R* SBL IN- 1 ) 

RETURN 

END 


2002 

2003 

2004 

2005 

2006 

2007 

2008 

2009 

2010 
2011 
20 12 
20 13 

2014 

2015 

2016 
2017 
20 18 
2019 

. 2020 
2021 
2022 
2023 
20 24 
20 25 
2026 
2027 



* 

SUBROUT INE TERPFIXI ,J, TA81, TAB2, TA83, TAB4,XI TAB ,EP) . 2028 

DIMENSION TABK 24),TAB2I 24) ,TAB3< 24) , TAB4I24) ,XITAB<24) 2029 

1 F I XI— .000 1 ) 2,2,10 2030 

2 GO TO I 3,4, 5,6),J 2031 

3 FP=2.53-2.439*AL0G< XI ) 2032 

GO TO 99 2033 

4 FP= 3. 54- 1 • 725*AL0G( »707l^XI ) 2034 

GO TO 99 2035 

5 FP=4.58-1.2195*ALDGl .5*XI ) 2036 

GO TO 99 2037 

6 FP= 10 . 1 2 2038 

GO TO 99 2039 

10 DO 12 N* l, 24 2040 

IFIXI-X ITABIN ) ) 11, 11, 12 2041 

11 NX=N 2042 

GO TO 13 2043 

12 CONT IN J E 2044 

13 T X= l X I-XlTAB(NX-l) )/l XITAB(NX)-XI TA8 I NX- 1 ) ) 2045 

TXl=l.-TX 2046 

GO TO ( 14, 15, 16, 17), J 2047 

14 FP«TXl*TABl(NX-lRrX*TABl{NXJ 2048 

GO TO 99 2049 

15 FP=TX 1 *T AB2I NX— 1 ) +TX*TAB 21 NX) 2050 

GO TO 99 2051 

16 FpaiTX 1*T AB3( NX— 1)*TX*TAB3(NX) 2 052 

GO TO 99 2053 

17 FP=TXl*TAB4rNX-mTX«TAB4fNX) 2054 

99 CONTINUE 2055 

RETURN 2056 

END 2057 
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SUBROUTINE SIMPINS, DX.ORD.F IND) . , 2058 

DIMENSION ORD( 50) '2059 

INTEGRATION OF NS ♦ 1 EQUALLY SPACED ORDINATE VALUES .^2060 

BY SIMPSON'S RULE. NS MUST BE EVEN • . <2061 

SUM = 0. 2062 

00 88 I =2. NS, 2 <: 2063 

SUM = SUM ♦ 2.*0RD( I-U ♦ A.*ORDII I , 206A 

FIND = DX *1 SUM - 0RDI1) ♦ ORD ( N S* 1 1 > / 3. , 2065 

RETURN , .... • . • , • 2066 

END . ? i- . 2067 
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* 

SUBROUT INE CORDXl N SBL ,NZ .ROBB ,SBL,X,XC> 2068 

2069 

BOUNOARY LAYER COORDINATES AND CORRESPONDING CHORDAL' 2070 

COORDINATES ARE COMPUTED HERE. 2071 

2072 

OIMENS ION SBL ( 300 ), XI 300) * XC I 300) 2073 

336 FORMAT! // 10X , 3 1 HITER AT ION TO COMPUTE XC FOR M -I5.32H DID NOT CONV 2074 

1ERGE IN 1000 STEPS. I 2075 

337 FORMAT! 1H1.2 5X , IHM, 20X , 1HS, 25 X, 1HX.24X.2HXC//) 2076 

338 FORMAT! 22X, 15. 3E25. 5) 2077 

M0UT=6 2078 

MX = NS BL ♦ NZ - 1 2079 

RZERO = RDBB/2. 2080 

XCI NZ ) =* -l. 2081 

DO 255 M=1,NZ 2082 

MM = NZ ♦ 1 - M 2083 

255 X ( M J = SBL (NZ I - SBL(MMI 2084 

DO 256 M=NZ,MX 2085 

MM = M *■ 1 - NZ 2086 

256 X { M ) = SBL l NZ ) * SBL I MM ) 2087 

DO 257 N=1,MX 2088 

IFINZ-M) 333,257,335 2089 

333 K * M > 1 - NZ 2090 

GO TO 334 2091 

335 K « NZ - M ♦ 1 2092 

334 XCIM) * -L. ♦ SBLIK) 2093 

l F I SBLI K ) — RZ'EKO T 341*341,342 2094 

341 XCIM) * -1. ♦ SBL(K)**2/(4.*R ZERO ) 2095 

342 CONTINJE 2096 

00 258 L»l. 1000 2097 

SAVE = XCIM) 2098 

CAL Cl « SORT I I l.+XC (M ) ) /RZERO ) 2099 

CALC2 = SORT ( T.-»( l . «-XCTM » I/RZERO) 2100 

XC(M) = XC(M)*CALCl*ISBL(X) - RZERO* (CALC 1*CALC2*AL0G (CALCl +CALC2 ) ) 2101 

D/CALC2 2102 

1 FI ABS I S AVE-XCI M I )- l.E-6) 257,257,258 2103 

258 CONTINJE 2104 

WRITE! MOUT , 336 ) M 2105 

257 CONTINJE 2106 

MRITE(MOUT,337) 2107 

DO 264 M*!1,MX 2108 

IFINZ-M) 261,261,262 2109 

262 K=NZ-M«-1 2110 

GO TO 263 2111 

261 K=H>I-N2 — 2112 

263 WRITEIMOUT, 338) M, SBLIK), XIM) , XCIM) 2113 

264 CONTINJE 2114 

RETURN 2115 

END 2116 
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SUBROUT INE PGRAOlMf X,UE,OXr,PRESS, SA, SB. SC , SR.SS) 

SUBROUTINE FOR CALCULATION OF PRESSURE GRADIENT AND 
DERIVATIVE COEFFICIENTS. 

DIMENSION X( 300), UE( 300.31' ' ' , '• . 

D1Z=X(M4-1I-X(M > • ' " ■' ’ 

02Z=X I M+2I-XIM ) ' y J' 

D21=X(*I*2I-X(M + 1) 

D1M1*XIM*1)-X(M-IJ ' 

DZM1=X l MJ— X(M— 1J 

XIM=DIZ/(D2Z*D21) 

ETAM= I./01Z-1./D21 L ' 

ZETAM=D2I/(D1Z*D2Z) \ 

PRESS »• (3.*UE(**l» II-A.*UEtMFiV2l+UEtH*I,3il/*2.*0XI)+UEIM+l ,1I»I 
IX IM*JE(M+2, I )+ETAM*UE4M*l,l )-ZETAM*U£ IM.II » 

S A= 1./01Z + 1./D1M1 . . . 

SBaDIMl/l D1Z*DZM1) 

SC=01Z/ID1M1*DZMI) 

SR=D1M 1/DZMl 
SS=01Z/DZMr' 

RETURN 


2117 
. 2118 

2119 

2120 
2121 
2122 
2123 
212A 

2125 

2126 

2127 

2128 

2129 

2130 

2131 

2132 

2133 

2 134 

2135 

2136 

2137 

2138 

2139 



* 

SU8R0UT INE TRANSIUPRIM, PRESS, THETA, REB .UC , NY, FLAM, XFLAM.LAMC) . 'V 2140 

C 1 ' ' ' - 2141 

C SUBROUTINE TO TEST FOR TRANSITION IN A LAMINAR BOUNDARY LAYER w \ 2142 

c • •; „ i; : ■ - - '’ 2143 

DIMENSION UC( 100,3), FLAMI 10), XFLAMCiO) 2144 

FIX) = .11746 - 1.0582E-3*X - 1 . I 022F-4* X* X , • : i _ . - V ' ■ 2145 

TKAY = PRESS*REB*THETA**2/UCINY ,2 ) .:•>-** r- 2146 

IFITKAY-.077) 2,2,99 • \>*)V ;.v’ 2147 

'.,2; IFIABSI TKAY)-. 0001) 3, 3,4 ' 2148 

3 ARG » IK AY* 72 • 48 ><.■•>? < .. - . . 2149 

GO TO 5 *, 2150 

4 ARG =0. 1 - 2151 

DO 6 N= 1, 1000 . . . ■ ■ • 2152 

SAVE = ARG . s ‘ j. j . :>■...?«■• 2153 

ARG = ARG - ( ARG* FI ARG )*,* 2- TK AY) / (F ( ARG ) *1.1 1,746— ARG*3. T746E— 3 — . A , 2154 

1RG*ARG*5. 5115E-4) ) « . . ,, *, -» .• . 2155 

IF( ABSI 1.-SAVE/ARGJ-1.E-6) 7, 7,6 ” ’ . 2156 

6 CONTINJE .. 2157 

7 I F I ARG* 11.) 8,8,5 : • 2158 

8 EF = 1.75 ' . „ . 2159 

GO TO 10 r • ' 2160 

5 • DO 15 N*l, 10 .» r : \ : 2161 

IF( ARG-XFLAMIN )) 24,24,15 2162 

24 NBAR = N 2163 

GO TO 16 2164 

15 CONTINJE 2165 

16 EF = FLAMI NBAR-1T*TARG-XFLAM( NBAR— 11 J*TFLAM (NBAR) -FLAMINBAR-l ) )/IX 2166 

1FLAMI NBAR )-XFLAM(NBAR- 1) ) 2167 

10 B = ,5*EF 2168 

A = 3.36*(UPRIM/UC(NY,2))**2 2169 

RTH » FI ARG)*( SORT! B*8*9860.*A)-B) /A 2170 

I FIREB*THETA-RTHI 99, 50, 50 2171 

50 LAMQ =0 2172 

99 CONTINJE 2173 

RETURN 2174 

END 2175 
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* 

SUBROUTINE C AP SI ITER , N , CAPO ,C AP H ,CAP J ,CAPK , SR ,SS ,SD ,SE ,SF , V I S C , V , U 2176 

1 1C ) 2177 

DIMENSION CAPGt 100),CAPH( 100) .CAPJtlOO) .CAPKIIOO) 2178 

DIMENSION VI SCI 100, 2) , VI 100,2 ),UC I 100,3) , SO (100) , SE 1 1 00 ) ,S F !1 00 ) 2179 

I F( ITER) 4,2,4 2180 

2 CAPG(N)= SR*V(N,1) - SS* V IN , 2 ) ; ' 2181 

CAPHI NI = SR*V ISCIN, U-SS*VISC(N,21 -- >’ ' * 2182 

CAPJI N ) = SR*( SOI N)*VISC(N + 1» 1 ) +SE IN)*VISC I N , 1 ) -SF ( N) *V I SC I N-l , 1 ) ) -S i 2183 

1S*I SDIN )*VISC( Nf 1, 2H-SEIN)*VI SC ( N , 2 )- SF I N) * VI SC IN-1 ,2) ) 2184 

C APK I N ) = SR*UC I N » 2)-SS*UC(N, 3) 2185 

GO TO 6 I 2186 

4 CAPG<N) = .5*(CAPG<M)+V(N, 1)1 2187 

CAPH(N) = .5*(CAPHm*VISC(N,H ) 2188 

CAPJI N ) = . 5*1 CAPJ(N)+SD(N)*VISCIN + 1,1)4SE I N> * VI SC ( N , 1 ) -SF I N) *VISCI N 2189 

1-1,1)) i 2190 

CAPK(N)=.5*ICAPKIN)4UC(N, 1) ) '< 2191 

6 CONTINJF 2192 

RETURN .’19) 

LNU ‘ 2 1V-. 

i 

4 


\ 
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* 

SUBROUI INE TERPI YIN.YBASE, VARY, NY, VALUE) • 2195 

2196 

SUBROUTINE FOR DETERMINING INTERPOLATED VALUE OF THE 2197 

FUNCTION VARY AT Y = YIN. 2198 

2199 

DIMENSION YBASEI 100) , VARY! 1001 2200 

1 FI Y IN-Y BASE! NY- 1 ) ) 2,3, 3 2201 

3 VALUE = V ARY ( NY ) 2202 

GO TO 10 2203 

2 DO 15 N- 1 , NY 2204 

IF( YIN-YBASEINI ) 24,24,15 2205 

24 N8AR=N 2206 

GO TO 16 2207 

15 CONTINUE . 2208 

16 D21=YBASEl NBAR J-YBASE ( NBAR-1 ) 2209 

D31*YBASE(NBAR*1»-YBASEINBAR-1) 2210 

032=031-021 2211 

D3A=YBASECNBAR*1 1-YIN 2212 

02 A=Y BASE! NBAR »- YIN 2213 

DA1=YIN-Y BAS El N BAR-11 2214 

V ALUE=D?A*D2A*VARYl NBAR-I 1/ (D21*D311*D3A*DA1 *VARYI NBAR1 / I D2i*D32 1- 2215 

102A*0A1*VARYI NBAR+l ) /( 031*032 1 2216 

10 CONTINUE 2217 

RETURN 2218 

END 2219 
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* 

SUBROUTINE YD I FF ( NY , ALPHA ,B ET A t GAMMA .DELTA , SD , SE , SF ,C2 ,C3,C4,Y) 2 2 20 

DIMENSION ALPHA! 100 ) # B ETA 1 1 00 ) ,GA MMA 1 1 00 ) .DELTAUOO) 2221 

DIMENSION SOI 100),SEI 1001.SFI 100) ,Y1100) 2222 

NV= NY -2 2223 

NVP1=NY«-1 2224 

DO 40 N = 2f NV ' 2225 

ALPHAm = 2«*!2»*Y(N)-Y!N-1)— YIN+ll) / !! Y! N+2 ) - Y I N— 1 ) ) * I Y I N+2 )— Y I N 2226 

lHll*(!Ut2)-Y(N))l 2227 

DELTA! N) =, 2.*! Y( N *2 ) +Yt N +1 J- 2.* Y IN) ) / ( ! Y( N«-2 > - Y ( N-l ) ) *1 Y < N*l ) -Y ! N" 22 28 

1-1) )*!YtN)-YtN-U>) 2229 

BET AIN) = (DELTACNJMYIN )-Y(N— 1 ) ) ** 3- ALPHA IN) *i Y IN*2 ) -YIN) ) **3 )/ !Y 2230 

1 I N+ 1 ) -T ( N ) ) **3 ! 2231 

GAMMAIN) = -ALPHAIN)-8ETAIN )- DELTA IN) ' 2232 

40 CONTINUE 2233 

00 39 N=2,NVPl 2234 

SDIN) = l YIN J-YIN-1) ) / II Y ( N + l )- Y I N— 1 ) )*IY!N+1)-Y!N) ) ) 2235 

SEIN) = 1 ./I YIN )-Y(N-l ) )- 1. / I Y'l N+D-YIN) ) 2236 

SPIN) = IYIN41)-YIN)>/IIY(N)-YIN-1»)*IYIN*1)-YIN-1))) 2237 

39 CONTINUE 2238 

C2 = Y( 3 )*Yl 4)/I Yl 2)*I Yl 3 » — YI 2) )* I YI4J-YI2) ) ) 2239 

C3 = -Y I 2 )*Yf 4')7"l Y'l 31* t Y( 4)-Y 1 3 ) I* ( Y( 31- Y(2) ) I 2240 

C4 = Y! 2)*Y! 3)/I YI 4)*I Y14)-YI 3) )*! Y!4)-YI2) ) ) 1 2241 

RETURN 2242 

END ' 2243 
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• . . * 

■ SUBROUTINE SETUP! LGO, M ,N V,REB ,X , Y ,UC , PRE SS ,GRAD ,DELT ,DI SP , THET A ,V D 2244 

ISC.MTRAN) 2245 

2 246 

SUBROUTINE FOR CALCULATION OF BOUNDARY LAYER THICKNESS, 2247 

DISPLACEMENT THICKNESS, MOMENTUM THICKNESS AND EDDY VISCOSITY. 2248 

. . - • 2249 

DIMENSION XI 300), Y( 100), UCIIOO, 3) , VI SC (100,2) .GRADIIOO) 2250 

RT R=SQRT I REB) , - , 2251 

NY = NV ♦ 2 2252 

U EDGE = ./J95*UCINY, 1) .. . . 2253 

DO 10 N= 1 , NV 2254 

(FIUEDGE-UClN+1, 1)) 41,41,10 . 2255 

41 NOELT = N 2256 

GO TO 20 2257 

10 CONTINJE , ■ . , 2258 

20 DELT = YINDELT )♦( UEDGE-UC (NDELT, 1 ) )* I YINDE LT41 ) -V I NDE LT ) ) / I UC I NDEL 2259 

IT 4-1,1 J-UCINDElT, l) ) j • i 2260 

SUM = 0. 2261 

DO 50 N*2, NY 2262 

50 SUM = SUM+I Y( N l-YI N- 1 ) )* I UC IN , 1 ) ♦UC.IN— 1 » 1 ) ) . . 2263 

DISP = (YINT)-.5»SUM/UC(NY,1) l/RTR 2264 

SUM * 0 . 2265 

U EDGE = UCINY.l) 2266 

DO 60 N*2 , NY 2267 

60 SUM = SUM*-(Y(N)-Y(N-1) )*( { UEDGE-UC IN,l))*UC(N,l)+( UEDGE-UC IN-1,1) ) 22 68 

1*UC( N-l , 1 ) ) 2269 

THETA =« .5*SUM/(RTR*UEDGE**2) 2270 

IFILGO) 53,53,56 2271 

53 NVP l=NV l 2272 

EASE = 1. 2273 

I FI M-MTRAN } 31,32,32 2274 

32 IFIMTRAN + 5-M) 31,31,33 2275 

33 EASE = IXTMT-XTWTRANl J/rXrMTRAN+SJ-XIMTRAN)) 2276 

31 CONTINJE 2277 

!NNER=0 2278 

FAC l = . 16*RTR*EASE 2279 

FAC2 = ,0168*UE0GE*DISP*REB*EASE 2280 

EFAC1 « -RTR/26. 2281 

EFAC2 « PRESS7RTR 2282 

TAUW = GRAD(1)/RTR 2283 

DO 160 N- 2* NVP l 2284 

ALTER = 1.+FAC2/I 1. *5. 5* < Y< N ) /DEL TJ**6) 2285 

IF! INNER) 402 "401, "402 2286 

402 V ISCIN, 1 ) * AL T ER 2267 

GO TO 160' 2288 

401 CONTINJE 2289 

TAUMY=TAUW-Y(N J*EFAC2 2290 

I FITAUMY ) 701,701,702 2291 

701 V ISCIN , U*1 2292 

GO TO 703 2293 

702 EX=YCN )*EFAC 1*SQR TITTORYT 2294 

VISCIN, l) = l.»FACi*Y(N)*Y(N)*ABS(GRAD(N))*(l.-EXP(EX))**2 2295 

703 IFIVISCINiD-ALTER) 160,160,521" 2296 

521 V ISCIN, l)*ALTER 2297 
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I NN ER= 1 2298 

160 CONT IN J E 2299 

S AV E= 1 • 2 300 

DD 162 N=2,NV 2301 

R AV E=V I SC ( N» 1 ) 2302 

V ISCIN, 1)=(V ISCIN*l, 1) +RAVE«-SAVE) /3. 2303 

162 S AVE=R AV E 2304 

56 CONTINUE 2305 

RETURN 2306 

END 2307 
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